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1  INTRODUCTION 


Steady  increases  in  large  scale  circuit  integration  indicate  that  the  Twenty-First 
Century  will  promise  significant  advances  in  High  Performance  Computing  (HPC) 
machinery.  Today,  one  may  obtain  desk-side  Linux  systems  containing  eight  processors 
(and  thirty-two  or  more  cores)  for  comparatively  reasonable  prices.  Moreover,  common 
laptop  systems  wield  significant  computing  power  with  central  processing  unit  (CPU) 
speeds  in  the  neighborhood  of  3.0  GHz  (maybe  more  by  the  time  this  report  is  certified) 
and  random  access  memory  (RAM)  storage  capability  in  hundreds  of  Gigabytes  (GB).  In 
the  realm  of  “Big  Iron”,  the  Department  of  Defense  (DoD)  High  Perfonnance  Computing 
(HPC)  Modernization  Office  recently  began  operating  clusters  each  with  tens  of 
thousands  of  cores,  and  the  Department  of  Energy  laboratory  community  has  even  larger 
systems.  These  developments  have  significant  implications  for  the  relatively  small 
Computational  Physics  research  community.  This  research  community  represented  by 
disciplines  such  as  high  energy  physics,  quantum  chemistry  and  computational  fluid 
dynamics  has  an  ever  increasing  need  for  computer  memory  and  for  parallel  processing 
speed. 


Computational  Fluid  Dynamics  (CFD)  has  drawn  on  HPC  resources  for  many 
years  to  help  with  aircraft  and  fluid  system  design.  Some  problems  like  high  Reynolds 
number  direct  numerical  simulations  are  still  computationally  inaccessible,  but  these 
situations  are  fewer  in  number  than  just  one  decade  ago.  For  instance,  we  routinely  solve 
problems  involving  the  large  eddy  simulation  (LES)  of  compressible  turbulence  with 
good  results.  Older  techniques  such  as  Reynolds-Averaged  Navier-Stokes  (RANS) 
simulation  now  teeter  on  the  brink  of  obsolescence.  Moreover,  massive  computing  power 
now  permits  us  to  invade  new  territory  previously  relegated  to  analytical  solutions 
supported  by  many  assumptions  and  highly  simplified,  under-resolved  computational 
studies.  Quantum  physics  now  benefits  widely  from  HPC  science  in  the  areas  of  quantum 
chemistry  and  molecular  dynamics.  These  areas  of  physics  now  impact  design 
engineering.  Although  it  occupies  only  a  very  small  part  of  the  research  community, 
detonation  physics,  a  close  relative  of  CFD,  can  benefit  handsomely  from  ever  more 
powerful  computational  techniques  and  equipment. 

1.0  Numerical  Detonation  Physics 

Numerical  Detonation  Physics  applies  many  of  the  same  computational 
techniques  employed  by  CFD.  The  primary  reason  is  because  detonations  are  powered  by 
the  propagation  of  the  detonation  wave,  a  powerful  shock  wave  that  transforms  the 
unreacted  explosive  into  detonation  product  species.  Like  the  shock  waves  encountered  in 
transonic  and  supersonic  flow,  detonation  waves  must  be  “captured”  in  the  material  field 
by  using  special  numerical  techniques.  Gas  phase  detonations,  e.g.,  the  explosive  bum  of 
acetylene  gas,  are  true  detonations  but  they  lack  some  of  the  complexity  associated  with 
the  detonation  of  condensed  (solid  or  liquid)  explosives.  Gas  phase  detonation  is  usually 
initiated  by  high  temperature.  It  follows  that  temperature  is  the  dominant  term  in  the 
reaction  rate  expression.  One  should  also  not  make  light  of  the  fact  that  we  actually  have 
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good,  quantitative  models  for  gas  phase  detonation  chemistry.  The  science  behind  the 
detonation  of  condensed  explosives  is  not  so  evolved. 

The  detonation  of  a  condensed  explosive  is  most  often  modeled  as  a  shock-driven 
process.  Macroscopic  observation  seems  to  indicate  that  a  shock  wave  is  often  required  to 
detonate  these  explosives.  Many  solid  explosives  simply  “burn”  when  exposed  to  a 
flame,  at  least  when  considered  over  relatively  short  time  periods.  Exposure  to  a  shock 
impulse  is  often  needed  to  initiate  the  run  to  detonation  for  an  explosive.  This  physics 
problem  is  complicated  greatly  because  of  the  smallness  of  scales  concerning  the 
detonation  wave.  The  detonation  wave  covers  a  thin  region,  a  fraction  of  a  millimeter  for 
most  ideal  or  Carbon-Hydrogen-Nitrogen-Oxygen  (CHNO)  explosives  like  Trinitro¬ 
toluene  (TNT).  The  head  of  the  detonation  wave  lies  at  the  entrance  to  the  detonation 
reaction  zone.  This  is  the  tiny  region  in  space  where  the  detonation  chemical  reactions 
take  place.  For  condensed  explosives,  we  do  not  know  these  chemical  reactions.  We 
know  only,  in  some  sense,  their  end  products,  and  if  we  detonate  two  like  samples  of  an 
explosive,  we  may  obtain  two  different  product  spectrums.  For  this  reason,  condensed 
explosives  are  relatively  crude  chemical  mixtures.  Still,  the  detonation  process  itself  may 
be  addressed  by  the  direct  application  of  the  conservation  laws  for  mass,  momentum  and 
energy.  This  same  approach  is  used  for  CFD  problems,  but  for  explosives  we  are  required 
to  apply  equations  of  state  for  both  the  unreacted  explosive  material  and  the  detonation 
products.  It  is  also  important  that  we  consider  heterogeneous  explosives.  These  materials 
contain  non-explosive  additives  like  plastic  binders  and  metal  particles.  In  future 
treatments  of  this  problem,  we  will  also  be  required  to  treat  the  material  behavior 
(material  strength  versus  applied  stress)  of  the  solid  explosive  in  response  to  shock 
excitation. 


1.1  A  Map  for  this  Report 

This  report  is  intended  to  assist  in  the  process  of  transitioning  detonation  physics 
algorithms  into  the  Farge  Eddy  Simulation  with  FInear  Eddy  Modeling  in  3  Dimensions 
(FESFIE3D)  multiphase  physics  computer  program.  The  discussions  that  follow  describe 
the  algorithms  applied  in  the  source  code  included  in  Appendix  A.  Although  these 
algorithms  are  tested  and  validated  to  some  extent,  it  is  nont  recommended  that  they  be 
coded  directly  into  FESFIE3D.  Rather,  the  Harten,  Fax  and  van  Feer  (HFF)  family  of 
algorithms  should  be  used  for  flux  difference  splitting  in  lieu  of  Roe’s  method.  Moreover, 
inhomogeneous  terms  in  the  equations  should  be  addressed  through  Strang  splitting. 1 

The  report  is  organized  as  follows.  In  Section  2,  we  describe  the  governing 
equations  for  the  detonation  problem  based  upon  the  work  of  Xu  et  al.2  Within  this  set  of 
equations,  we  add  the  terms  coupling  the  detonation  flow  field  to  the  particle  field.  We 
show  that  reaction  rate,  particle  coupling  and  geometric  effects  may  be  incorporated  as 
source  terms.  The  equations  of  state  used  for  the  solid  explosive  and  for  the  detonation 
products  are  also  presented  in  this  section.  The  advective  terms,  of  critical  importance  in 
the  shock-capturing  scheme,  are  clearly  delineated.  Section  3  describes  the  eigen- 
structure  for  the  system  of  governing  equations.  The  flux  Jacobian  matrix  is  developed 
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for  the  reactive  Euler  equations  adapted  for  a  real  gas  equation  of  state.  Then  we  develop 
a  set  of  eigenvalues  and  eigenvectors  needed  in  order  to  accurately  capture  the  detonation 
wave.  In  Section  4,  we  discuss  the  overall  numerical  scheme  and  temporal  discretization 
procedure  used  in  our  detonation  computer  program.  We  also  discuss  the  development  of 
the  numerical  flux  vector  in  detail.  Section  5  contains  the  terms  governing  the  motion  of 
Lagrangian  particles  including  the  drag  laws.  In  Section  6,  we  provide  the  results  for 
three  example  calculations.  After  perfonning  a  calculation  to  verify  proper  code 
perfonnance,  we  simulate  the  detonation  of  a  spherical  mass  of  HMX  loaded  with  metal 
particles.  We  show  a  series  of  detonation  waveforms  for  this  explosive,  and  we  go  on  to 
include  the  resulting  particle  trajectories  and  velocities.  We  also  make  some  basic 
comparisons  between  the  results  produced  by  our  computer  program  to  archival 
explosive  performance  data  for  HMX.  Finally,  in  Section  7,  we  draw  several  important 
conclusions  from  our  development.  We  also  make  recommendations  for  follow-on  work 
needed  to  support  the  installation  of  detonation  physics  algorithms  in  LESLIE3D. 


Distribution  A.  Approved  for  public  release,  distribution  unlimited.  (96ABW-201 1-0548) 

3 


2  GOVERNING  EQUATIONS 


To  address  the  detonation  problem,  we  follow  a  body  of  research  documented  in 
the  general  scientific  literature.2  By  doing  so,  we  can  escape  some  of  the  uncertainties 
associated  with  the  older  programmed  burn  detonation  models.  We  do  make  a  departure 
from  the  core  reference  in  that  our  development  disregards  the  issue  of  compaction  in  the 
solid  explosive.'  Instead,  it  is  assumed  that  our  explosive  is  a  solid  mass  at  or  near  the 
theoretical  maximum  density.  The  present  approach  allows  the  reaction  zone  to  be  clearly 
resolved  within  the  limitations  of  the  grid  refinement.  As  a  result,  the  forces  applied  to 
particles  may  be  resolved  more  accurately. 

2.1  The  Reactive  Euler  Equations 

The  reactive  Euler  equations  are  frequently  used  to  represent  detonation  flow 
fields  based  upon  a  reaction  progress  equation  and  a  mixture  equation  of  state."  The 
equations  for  the  conservation  of  mass,  momentum,  energy  and  reaction  progress  may  be 
readily  expressed  in  vector  form.  The  equation  for  a  detonation  field  set  in  one  space 
dimension  may  be  written  as 


—  +  — =  SG+Sfe+SP  (2.1.1) 

dt  8x  G 

where 

U  =  [p,  pu,  E,  pAf  (2.1.2) 

is  the  vector  of  conserved  variables,  and 

F  =  \pu,  pu2  +  P,  u(E  +  P),puAf 

(2.1.3) 

is  the  flux  vector.  Also, 


SG  =  -—\pu,  pu2,  u(E  +  P),  puA ]7 

(2.1.4) 

X 

[0,0,0  ,pr]T 

(2.1.5) 

SP=[0,Fs,Qs,0f 

(2.1.6) 

We  may  also  write  the  total  energy  per  unit  volume  as 

E  =  pe  +  ^u2  (2.1.7) 
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where  e  is  the  internal  energy  per  unit  mass.  The  equation  of  state  may  be  written  in  the 
general  fonn 


P  =  P{p,  e.  A) 


(2.1.8) 


where  A  is  the  reaction  progress  variable. 

Vectors  SG,  and  Sp  contain  source  tenns;  as  we  have  shown,  these 
nonhomogenous  tenns  are  kept  on  the  right  hand  side  of  the  reactive  Euler  equations  and 
may  be  treated  independently  from  the  advective  terms.  Vector  Sc  contains  the 
geometric  source  terms  that  allow  the  system  to  be  configured  for  planar,  cylindrical  or 
spherical  one-dimensional  flow.  To  adapt  (2.1.1)  for  planar  flow,  we  need  only  set  /  =  0 
in  (2.1.4).  We  may  adapt  (2.1.1)  for  cylindrical  or  spherical  one-dimensional  flow  by 
setting  ;  =  1  or  /' =  2,  respectively.  Vector  S  contains  the  reaction  rate  source  tenn 

governing  the  rate  of  progress  for  the  detonation  reaction.  The  reaction  rate  r  may  be 
written  in  many  different  forms  depending  on  the  explosive.4  The  tenn  we  have  chosen  to 
use  for  HMX  may  be  written  as 


r  =  k 


f  \N 

'  P  ' 
KPcjJ 


(1-^r 


(2.1.9) 


where  PC1  is  the  Chapman- Jouquet  pressure  for  HMX;  k ,  N  and  v  are  constants 
chosen  to  fit  experimental  data.5  Note  that  this  reaction  rate  law  is  dependent  upon  both 
pressure  and  reaction  progress.  The  source  term  vector  SP  has  been  added  to  the  system 
by  the  author.  It  represents  the  dynamic  coupling  between  the  detonation  products  and  a 
field  of  discrete,  massive  Lagrangian  particles.  The  coupling  is  based  upon  both 
momentum  and  thennal  effects.6  The  specific  forms  of  the  coupling  terms  are  presented 
in  a  later  section. 

2.2  Mixture  Equations  of  State 

For  the  detonation  problem,  relevant  equations  of  state  are  cast  in  the  form  of 
(2.1.8).  This  fonn  is  complicated  since  pressure  varies  as  a  function  of  density,  internal 
energy  per  unit  mass  and  reaction  progress.  In  this  analysis,  the  reaction  progress  variable 
is  analogous  to  a  species  mass  fraction  commonly  used  in  reacting  gas  flows.  Moreover, 
it  is  used  to  compute  the  specific  internal  energy  for  the  detonating  mixture  by  fonning  a 
weighted  sum  of  the  equation  of  state  (EOS)  for  the  solid  explosive  and  the  EOS  for  the 
detonation  products.  The  resulting  expression  for  specific  internal  energy  is  called  the 
mixture  EOS.2  Our  governing  equations  (2.1.1),  discretized  in  accordance  with  the  finite 
volume  method,  rely  upon  the  mixed  cell  approach.  Each  flow  cell  is  assumed  to  contain 
a  mixture  -  part  solid  explosive  and  part  detonation  products.  The  mixture  fraction  is 
given  by  the  reaction  progress  variable  A  ,  and  A  is  defined  as  the  mass  fraction  of  the 
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detonation  products  in  the  cell.  The  density  within  a  cell  is  the  sum  of  the  densities  for  the 
solid  ( 5  )  and  gas  ( g  )  phases,  respectively,  i.e., 


P  =  Ps+Pg 


so  A  is  given  by 


and 


^  =  1  -A 
P 


(2.2.1) 


(2.2.2) 


(2.2.3) 


Hence,  we  have  that  A  is  the  mass  fraction  of  the  gas  (detonation  products)  phase.  We 
also  assert  that  the  internal  energy  for  a  given  finite  volume  cell  may  be  expressed  as 

e  =  Aeg+(1-A)es  (2.2.5) 


where  eg  and  es  are  the  specific  internal  energies  for  the  gas  and  solid  phases, 

respectively.  This  mixing  rule  differs  from  the  archived  approach  based  upon  specific 
volume,  but  to  date,  we  have  not  been  successful  in  applying  Xu’s  closure.7  Assume  the 
same  pressure  for  both  phases  with  each  phase  having  its  own  equation  of  state,  i.e., 

eg=eg{Pg,P)  (2.2.6) 

es=es(ps,P)  (2.2.7) 

with  pg  and  ps  given  by  (2.2.2)  and  (2.2.3). 


2.3  Solid  Explosive  Equations  of  State 

In  the  previous  section,  we  showed  that  one  part  of  our  mixture  EOS  represents 
the  solid  explosive.  In  the  discussions  that  follow,  we  apply  two  different  fonns  of  an 

o 

EOS  originally  developed  by  Hayes.  The  first  fonn  of  this  EOS  (Hayes-I)  works  very 

2 

well  for  mechanical  effects.  The  Hayes-I  EOS  is  given  as 


P-P 

es(ps,P)  = - 51 

8 


( 

h 

V 


r>  h 

f  \ 

f  \ 

N- 1 

r  \ 

A 

o| 

Q.I 

1 

+  t4 ‘ 

Ps 

~(N~  1) 

o  1 

Oj 

Q.I 

1 

-1 

Ps  0  j 

^  Ps  ) 

\  Ps  0  ) 

V  Ps  ) 

(2.3.1) 


where 
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8  Tq  ps  q 

(2.3.2) 

t  _  Cvs  ro  8 

Pso 

(2.3.3) 

t  -  H 1 

4  psQ  N(N-T) 

(2.3.4) 

In  equations  (2.3.1)  through  (2.3.4),  P0,  T0  and  ps0  are  the  ambient  pressure, 
temperature  and  unloaded  solid  density.  ro  is  the  Gruneisen  parameter,  and  Cvs  is  the 

constant  volume  specific  heat  for  the  solid.  Hl  and  N  are  parameters  used  to  fit  the  EOS 

2 

to  data.  Table  1  lists  all  of  the  required  parameters  for  this  EOS. 


Table  1  -  Hayes  EOS  Data  for  HMX 


Hi 

1.3  x  1010  N/m2 

N 

9.8 

Cvs 

1.5  x  103  J/(KgK) 

To 

1.105 

Po 

101325  Pa 

PsO 

1.9  x  103  Kg/m3 

To 

300  K 

The  second  form  of  the  Hayes  EOS  (Hayes-II)  functions  well  mechanically  but 
also  incorporates  temperature.  The  Hayes-II  EOS  is  given  as 


es(Ps>p)  =  - 
8 


P-P0~  —  < 

0  N 


f  \N 

Ps 

Ps  0  y 


-1 


f 

h 

V 


Pso  A 


y 

PsO 

Ps  y 


f  \ 

N- 1 

(  \ 

' 

Ps 

-(NY) 

|  \  PsQ 

-1 

V  PsO  ) 

{  Ps  J 

(2.3.5) 


This  version  of  the  Hayes  EOS  may  be  derived  by  using  Reference  1;  however, 
additional  terms  are  incorporated  in  (2.3.5)  to  match  the  behavior  of  (2.3.1)  at  ambient 
pressure.  The  temperature  of  the  solid  explosive  is  given  by 


T(p„P)  =  - 


-  H'  ■ 
0  N 


v 


f  \N 

_Ps_ 

PsO  y 


+  Tn 


(2.3.6) 


Distribution  A.  Approved  for  public  release,  distribution  unlimited.  (96ABW-20 1 1  -0548) 

7 


Together,  equations  (2.3.5)  and  (2.3.6)  constitute  a  complete  equation  of  state  for  a  solid 
explosive.9  These  equations  use  the  same  data  as  is  listed  in  Table  1  for  HMX.  The 
Hayes-II  EOS  also  perfonns  very  well  in  one-dimensional  detonation  studies  for  solid 
HMX. 


2.4  Detonation  Products  Equation  of  State 

As  equation  (2.2.5)  indicates,  part  of  the  mixture  EOS  must  address  the  gaseous 
products  resulting  from  the  detonation  of  the  solid  explosive.  For  the  purposes  of  this 
work,  we  have  selected  the  .Tones-Wilkins-T.ee  (TWL)  EOS.1  The  TWL  EOS  is  somewhat 
controversial,  but  nevertheless,  it  is  widely  applied  in  hydrocodes.  Also,  many  explosives 
have  been  characterized  for  this  EOS.  We  apply  the  TWL  EOS  in  the  following  fonn. 


e,(p,,P)  = 


cop. 


P-A 


c°Pg 

A 


exp 


PSJ 


-B 


cop. 


R 


exp 


2  J 


R 2 
P* 


-<2  +  e0  (2.4.1) 


where  A ,  B  ,  a> ,  Rx  and  R2  are  coefficients  produced  by  curve-fitting  for  the  explosive 
under  consideration.  Also,  note  that 


A  =  *1  Pso  .  (2-4.2) 

and 

A  =  Ri  Pso  •  (2-4-3) 

Q  is  the  heat  of  detonation  for  the  explosive,  and  e0  is  the  reference  value  for  specific 
internal  energy.  There  is  no  firm  rule  for  determining  e0 ,  but  we  will  define  e0  as 

e0  =  Cvg  T0 .  (2.4.4) 


Table  2  -  TWL  Coefficients  for  HMX 


Ri 

4.2 

r2 

1.0 

CO 

0.3 

A 

7.783  x  1011  Pa 

B 

7.071  x  1010  Pa 

Cvg 

(1.1  -  0.28x1  O'3  ps0)  x  103  T/(Kg  K) 

Q 

[7.91  -  4.33  (10-3ps0  -1.3)"  -  0.934  (10-3ps()  -1.3)] 
x  106  T 

C  is  the  constant  volume  specific  heat  for  the  detonation  products.  The  data  used  for 

9 

HMX  in  the  TWL  EOS  is  listed  in  Table  2.  For  the  studies  performed  later  in  this  work, 
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mixture  EOS. 
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3  SYSTEM  EIGEN-STRUCTURE 


3.1  Flux  Jacobian  Matrices 


Capturing  the  structure  of  the  detonation  wave  constitutes  a  difficult  numerical 
issue  involving  the  discretization  of  the  advective  tenn  ;  where 


t  0F 
A  = - = 

au 


a f, 

dF ; 

dF, 

dF ; 

dp 

d(pu) 

dE 

dx 

8F2 

dF2 

dF2 

dE2 

dp 

d(pu) 

dE 

dx 

dF3 

dF3 

dF3 

dF3 

dX 

d(pu) 

dE 

dX 

dF4 

dFA 

dF4 

dF4 

ds l 

d(pu) 

dE 

dx 

(3.1-1) 


is  called  the  flux  Jacobian  matrix.  The  tenn  Ft  simply  denotes  the  ith  element  of  the  flux 

vector  F.  Equation  (3.1.1)  is  already  annotated  with  the  specific  elements  of  U .  It  is 
important  to  note  that  our  equation  of  state  is  cast  in  a  general  fonn,  so  the  calculation  of 
the  specific  elements  of  (3.1.1)  is  made  more  complicated.  The  method  for  calculating 
these  matrix  entries  relies  heavily  on  the  derivatives  of  pressure  taken  with  respect  to  the 
conservative  variables.10  For  convenience,  the  pressure  derivatives  for  this  Jacobian  are 
given  below.  For  the  three-dimensional  case,  the  detailed  derivation  of  these  pressure 
derivatives  is  presented  in  Reference  1 1 .  For  pressure  given  in  the  form  of  (2. 1 .8),  let 


pp  = 


dP) 


;p  = 


ap 

yfejpjL 


;p,= 


'dP' 


p,e 


then  we  may  write  the  pressure  derivatives  as 


f  dp^ 


dP)Pu 


=pp+pe 


f  2 

U 


\^yjPu,E,P* 


p  P  ) 


p 


dP 

Kd{pu))p^p, 


=  ~U~Pe 

P 


(3.1.2) 


(3.1.3) 


(3.1.4) 


'  6P_\  _PL 
)  P,Pu,pX  P 


(3.1.5) 
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J3  =  (H -u2)^  +  A^-  (3.1.9) 

P  P 

and  the  frozen  speed  of  sound,  a  ,  is  given  by 

P  P 

a2=Pp+ — T-  (3-1-10) 

P 

The  derivation  for  this  speed  of  sound  is  also  archived. 1 1 

We  can  also  define  a  vector  of  non-conservative  variables  for  the  reactive  Euler 
equations  as  V  ,  where 

y  =  [p,u,P,Af.  (3.1.11) 

As  you  may  surmise,  the  governing  equations  may  also  be  written  in  terms  of  the  non¬ 
conservative  variables,  and  we  may  define  a  non-conservative  flux  Jacobian  matrix  A 
such  that1 1 
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(3.1.12) 


u 


A  = 


0 

0 

0 


p 

u 

pa 2 
0 


0  0 

0 

u  0 

0  u 


The  derivation  of  the  non-conservative  reaction  progress  is  a  simple  exercise.  Observe 
that  the  conservative  form  of  this  equation  is  written  as 


d(pA)jd(pA) 
dt  dx 


(3.1.13) 


We  may  expand  (3.1.13)  as  follows. 


A 


dp  |  d{pu ) 
v  dr  dx 


^  dA  dA 
+  P-  +  PU- 


dt 


dx 


pr 


(3.1.14) 


The  first  tenn  in  (3.1.14)  vanishes  since  it  is  just  a  scalar  multiple  of  the  continuity 
equation  (component  one  of  2. 1 . 1),  so  we  obtain 


dA  dA 

- h  u  —  =r 

dt  dx 


(3.1.15) 


as  the  non-conservative  reaction  progress  equation. 

3.2  Eigenvalues 

The  eigenvalues  of  the  flux  Jacobian  matrix  contain  important  infonnation  on  the 
physics  of  our  detonation  problem.  We  think  of  any  fluid  mechanics  problem  (as  well  as 
most  solid  mechanics  problems)  in  terms  of  interacting  waves.  The  detonation  problem 
can  be  decomposed  into  a  set  of  characteristic  waves.2  The  speeds  at  which  these  waves 
propagate  are  given  by  the  eigenvalues  of  the  flux  Jacobian  matrix.  For  any  square 
matrix  A  ,  the  eigenvalues  are  defined  as  the  set  of  numbers  C,  such  that 


A-£I =0 


(3.2.1) 


where  1  is  the  identity  matrix.  We  may  note  that  the  conservative  matrix  (3.1.7)  is 
heavily  populated,  so  it  is  very  difficult  to  obtain  the  eigenvalues  by  using  (3.2.1). 
Fortunately,  the  non-conservative  matrix  (3.1.12)  is  a  simpler  fonn  mathematically 
equivalent  to  (3.1.7),  so  these  matrices  must  have  the  same  eigenvalues.11  Using  (2.3.1), 
the  eigenvalues  of  (3. 1 . 12)  are  easily  shown  to  be 

e  {u -a,  u,  u,  u  +a}  (3.2.2) 
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Note  that  u  is  an  eigenvalue  of  multiplicity  two,  so  there  are  two  waves  with  speed  u  , 
i.e.,  the  entropy  and  reaction  progress  waves  both  propagating  at  the  flow  velocity.  The 
remaining  two  distinct  eigenvalues  C,  =  u  ±  a  denote  acoustic  waves.  The  dynamics  of 
the  detonation  process  may  be  described  through  the  interactions  of  characteristic  waves, 
but  to  completely  describe  these  waves,  we  must  detennine  the  eigenvectors  for  the 
detonation  problem. 

3.3  Eigenvectors 

In  order  to  detennine  the  characteristic  waves  for  (2.1.1),  we  must  detennine  the 
eigenvectors  for  the  conservative  Jacobian  matrix  (3.1.7).  When  we  use  the  tenn 
eigenvector,  in  this  case,  we  are  referring  to  a  right  eigenvector .10 

Definition'.  Given  a  matrix  A  e  C (n  x  n)  with  a  set  of  eigenvalues  e  C ,  i  =  1, . . . , n ,  we 
define  the  right  eigenvector  r;  e  C(n)  associated  to  the  eigenvalue  C!  l  such  that 

Ar,=^.r,  (3.3.1) 

Equation  (3.3.1)  is  useful  in  that  it  tells  us  how  to  find  right  eigenvectors.  To  find  a  right 
eigenvector  for  (3.1.7)  associated  to  an  eigenvalue  C,  ,  we  first  define  the  components  of 
right  eigenvector  r  .  Let 

r  =  (ovo2,o3,u4)t  (3.3.2) 

Now  we  apply  (3.1.7)  and  (3.3.1)  to  create  a  linear  system  of  equations  in  the 
components  of  r . 

0 

a 2  -  u 2  -  / 3 

u(a2  -  H  -  J3) 

-u  A 

The  system  (3.3.3)  directly  leads  to  a  system  of  four  eigenvector  equations.  The 
eigenvector  equations  do  not  have  a  unique  solution;  in  fact,  they  have  an  infinite  number 
of  solutions,  so  care  is  required  in  structuring  prospective  choices  for  the  components  of 
r  to  design  a  proper  numerical  treatment  for  the  problem.  Also,  it  is  important  to  observe 
that  the  number  of  linearly  independent  eigenvectors  must  be  same  as  the  order  of  the 
system.  For  this  detonation  problem,  the  Jacobian  matrix  is  of  the  fourth  order,  so  we 
must  detennine  four  linearly  independent  eigenvectors  even  though  we  have  only  three 
distinct  eigenvalues;  the  eigenvalue  u  is  repeated. 

We  begin  the  process  of  determining  some  specific  eigenvector  components  by 
extracting  the  first  eigenvector  equation  from  (3.3.3),  i.e., 


2-  — 

PJ 


H  -  —  P, 


P 

A 


0 

K 

P 

1  +  ^L 

V  P) 

0 


0 

P± 

P 

-P 

1  A 

P 


p 

hi 

o2 

II 

C\J 

3 

v3 

o4 

_°4_ 

(3.3.3) 


Distribution  A.  Approved  for  public  release,  distribution  unlimited.  (96ABW-20 1 1  -0548) 

13 


»2=C»1 


(3.3.4) 


We  may  satisfy  equation  (3.3.4)  by  choosing 

=1;  v2=^  (3.3.5) 

Equation  (3.3.5)  may  be  used  in  (3.3.3)  to  produce  the  remaining  three  eigenvector 
equations 

a2-u2-j3  +  [  2 u-C~-Pe  W  +— y3+  —  =  0  (3.3.6) 

l  P  )  P  P 

(  2  \  (  \ 

u(a2  -  H  -  j3)  +  £  H-—  +  i-C  +  -Pe  u3+-PAv  4=0  (3.3.7) 

v  P )  v  P  )  P 

-uA  +  £  A  +  (u-£)v4=0  (3.3.8) 

Based  upon  (3.3.5),  we  may  produce  the  eigenvector  associated  to  eigenvalue  £  =  u  .  Set 
C,  =  u  in  (3.3.8),  and  we  see  that  this  equation  is  trivially  satisfied  with  no  restrictions  on 
v4 .  Now  we  set  £  =  u  in  (3.3.7)  and  (3.3.8);  by  simplifying,  we  can  show  that  both  of 
these  equations  reduce  to  the  same  equation,  i.e., 

u2  P  P 

a2  -j3-—Pe  +  ^o,+^u4=  0  (3.3.9) 

P  P  P 

Since  there  are  no  restrictions  on  o4  ,  we  may  freely  choose  v4  and  solve  for  v3 . 

u3=H-^  +  ^-(A-u4).  (3.3.10) 

1  e  1  e 

By  cleverly  choosing  the  value  of  o4  ,  we  produce  two  linearly  independent  eigenvectors 
associated  to  the  eigenvalue  £  =  u  .  If  we  set  v4  =  0 ,  we  obtain  the  eigenvector 

(  n  2  P  Y 

r  =  \,u,H  ~^  +  ^A,0  (3.3.11) 

pp 

V  re  rX  7 

Alternatively,  we  obtain  a  second  eigenvector  by  setting  o4  =  1 ,  so 

(  2  p  y 

r  =  +  (3.3.12) 

V  Pe 
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We  may  also  obtain  the  eigenvector  associated  to  eigenvalue  £  =  u  +  a;  by 
returning  to  equation  (3.3.4),  let  us  choose 

l>j=1  ;  u2=u  +  a  (3.3.13) 

By  substituting  (3.3.13)  into  (3.3.8),  we  may  show  that 

v4=A  (3.3.14) 

We  can  produce  another  eigenvector  equation  associated  with  this  eigenvalue  by  using 
(3.3.14)  and  setting  £  =  u  +  a  in  (3.3.6).  By  doing  so  and  solving  for  t>3  ,  we  have  that 

t>3  =H+ua  (3.3.15) 

One  may  show  that  (3.3.13),  (3.3.14)  and  (3.3.15)  satisfy  (3.3.7),  and  the  eigenvector 
associated  to  eigenvalue  C,  =  u  +  a  is 

r  =  (1  ,u  +  a,  H +  ua,  A)T  (3.3.16) 

We  may  derive  the  eigenvector  associated  to  eigenvalue  C,  =u  —  a  by  the  same 
procedure.  We  consider  (3.3.4)  and  then  set 

=1;  v2  =u-a  (3.3.17) 

Equation  (3.3.8)  can  be  applied  to  again  obtain  the  result  (3.3.14).  By  substituting 
(3.3.17)  and  (3.3.14)  into  (3.3.6),  we  can  solve  for  u3  ,i.e., 


u3=H-ua.  (3.3.18) 

Subsequently,  one  can  show  that  (3.3.17),  (3.3.18)  and  (3.3.14)  satisfy  equation  (3.3.7). 
Hence,  the  eigenvector  associated  to  eigenvalue  £  =  u  —  a,  may  be  written  as 


r  =  (1,  u -a,  H -ua,  A)T 


(3.3.19) 


Equations  (3.3.11),  (3.3.12),  (3.3.18)  and  (3.3.19)  are  the  eigenvectors  for  the  reactive 
Euler  equations  in  one  dimension.  We  can  form  R  ,  the  matrix  of  right  eigenvectors,  by 
allowing  each  eigenvector  to  form  a  column  of  this  matrix.  Hence, 


1 


1 


R 


u  -  a  u  u  u  +  a 

2  2 

H-ua  +  H-^-  +  ^(A- 1)  H+ua 

P  P  PP 


A 


0 


1 


A 


(3.3.20) 
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It  is  a  straightforward  although  tedious  exercise  to  show  that  R  ,  the  determinant  of  R , 
is 


(3.3.21) 


So  far,  our  development  of  the  eigen-structure  for  the  reactive  Euler  equations  closely 
coincides  with  Glaister’s  derivation  perfonned  for  the  real  gas  equation  of  state.10  From 
(3.3.21),  we  can  see  that  our  eigenvectors  are  well-defined  and  constitute  a  non-singular 
system  for  realistic  values  of  density  and  the  speed  of  sound  with  Pe  ^  0 .  As  a  result,  R 
is  invertible  under  the  same  conditions,  and  we  can  calculate  the  matrix  of  left 
eigenvectors  L  with  L  =  R  ,  and  by  using  the  adjoint  matrix  for  R  (the  transpose  of 
the  matrix  of  cofactors)  in  conjunction  with  the  definition  of  the  inverse  matrix,  we  have 
that 


H  -u2  - (u  +  a)  +  —  / i 


le  J 


2a 


(1  -X)(u2 -H)-j{pa2 -P^-V>) 


\ 

a 

/9fl' 
U  +  — — 

) 

'j 

l  Pe) 

R 

f 

2a  A 

V 

a 


uz-H  +  y(pa2-APA) 


H  -u2  +  ^-(u-a)  +  —  A 

Pe  Pe  J 


2 u  a(A  - 1) 
-  2  u  a  A 


\ 

a 

’  pa ^ 

u - 

) 

l  Pe) 

-a 


-a- 


2u(1  ~  A)  2a 
2 a  A  2 a\ 


-a 


V  P.  P, 

'  e'+'Lx 

V  Pe  Pe  J 


-a 


J 


(3.3.22) 


Each  row  of  the  matrix  shown  in  (3.3.22)  is  a  left  eigenvector  for  the  Jacobian  matrix 
found  in  (3.1.7). 

Although  we  have  not  yet  presented  explicit  fonns  for  the  pressure  derivatives, 
we  have  accomplished  a  great  deal  of  work  in  this  section.  Equations  (3.2.2),  (3.3.20)  and 
(3.3.22)  offer  a  complete  description  of  the  structure  of  the  eigen-space  associated  with 
the  flux  Jacobian  matrix  A  shown  in  3.1.7.  Moreover,  we  can  fonnulate  a  special 
similarity  transformation,  i.e., 
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A  =  RAL 


(3.3.23) 


or 


A  =  LAR 


(3.3.24) 


and 


u  -a 


0  0  0 

u  0  0 

0  u  0 

0  0  u  +  a 


(3.3.25) 


is  the  diagonal  matrix  of  eigenvalues.11  Recall  that  matrix  L  is  the  inverse  of  R.  Our 
discussion  of  the  numerical  physics  behind  Roe’s  scheme  for  the  reactive  Euler  equations 
is  now  complete.  The  Roe  formulation  is  quite  important  from  the  theoretical  standpoint, 
but  this  method  is  difficult  to  implement  for  two  or  more  non-Cartesian  space 
dimensions.  Fortunately,  other  flux -based  discretization  methods  such  as  the  Harten,  Lax 
and  van  Leer  (HLL)  family  of  schemes  can  easily  be  applied  to  this  problem.  Moreover, 
these  methods  do  not  require  the  calculation  of  pressure  derivatives  (yet  to  be  discussed) 
for  the  mixture  equation  of  state.  This  fact  affords  greater  of  ease  of  calculation  for  a 
production  numerical  scheme. 
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4  BUILDING  THE  NUMERICAL  SCHEME 


In  this  section,  we  pull  together  all  of  the  aspects  of  detonation  physics  and 
mathematics  discussed  in  preceding  sections  and  dedicate  our  efforts  to  the  solution  of 
our  benchmark  problem  -  simulating  the  detonation  of  a  finite  sphere  of  HMX.  In  order 
to  accomplish  this  goal,  we  begin  by  presenting  detailed  pressure  derivatives  for  our 
mixture  equation  of  state.  Then  we  discuss  the  details  associated  with  our  chosen 
numerical  integration  scheme  including  formulation  of  the  numerical  flux  vector. 

4,1  Pressure  Derivatives 

The  purpose  of  this  subsection  is  to  document  formulas  for  the  pressure 
derivatives  (3.1.2)  of  the  mixture  equations  of  state.  These  derivatives  must  be  computed 
under  the  support  defined  by  the  set  of  primitive  variables.11  In  this  work,  we  consider 
two  mixture  equations  of  state.  The  first  mixture  EOS,  called  the  Hayes-I/JWL  EOS  is 
given  by  substituting  (2.3.1)  and  (2.4.1)  into  (2.2.5).  The  second  mixture  EOS,  referred 
to  as  the  Hayes-II/JWL  EOS,  is  created  by  substituting  (2.3.5)  and  (2.4.1)  into  (2.2.5). 
Either  mixture  EOS  consists  of  a  lengthy  formula,  so  to  promote  brevity  in 
documentation,  we  can  relate  the  two  mixtures  equations  of  state  to  one  another.  If  we 
look  carefully  at  the  Hayes-I  and  Hayes-II  formulas,  (2.3.1)  and  (2.3.5),  respectively,  we 
see  that 


a  i  H  \ 
e,  =  e 


s  s 


gN 


'p_'N 

v  A)  y 


-1 


(4.1.1) 


These  expressions  for  the  internal  energy  of  the  solid  explosive  differ  by  only  one  term. 
The  Hayes-I/JWL  mixture  EOS  may  be  written  as 


4  =0-^)  el  +^e 


(4.1.2) 


Hence,  by  using  (4. 1 . 1),  we  may  write  the  Hayes-II/JWL  mixture  EOS  as 


4=(i-2)e;-Mzy) 

gN 


Po  . 


+  A  6  „ 


(4.1.3) 


where  we  have  used  (2.2.3).  A  general  formula  for  the  Hayes- X/JWL  mixture  EOS  may 
be  written  as 


PK  -  c\-2  V7  -fiK  —  ^  ^ 

eM  0  -4)e  djj 

gN 


/  \  N 

Po  . 


1 


+  A,  e  n 


(4.1.4) 


Accordingly,  equations  (2.3.1)  through  (2.3.4)  may  be  used  to  expand  (4.1.4)  and  obtain 
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eM  =  P  D-  (3 


-A 


where 


o| 

1 

1 

+  tA  0  “ 

-xf 

(JL 

V  P ) 

{  Po) 

AM 


~hi  1-A) 


f  1  2^ 

r  r,  ^ 

f  1  2^ 

f  r2 

- — 

exp 

-B 

- — 

exp 

— 

ycop  r^j 

l  xpj 

{cop  R2J 

1  xp) 

~(Q  +  e0)A  (4.1.5) 


-8 


kH^-X) 


gN 


f  /j  0,  \N 

0  -X)p 

Po  / 


-1 


„  1-2  1 
D  = - +  — 


g  cop 

p 

e=t 3 — 1 

Po 

/?  =  6?+(A-1)r4 


1 5  =?4+  — 
<? 


(4.1.6) 

(4.1.7) 

(4.1.8) 

(4.1.9) 


Equation  (4. 1 .5)  may  be  solved  for  pressure,  i.e., 


P  =  - 

D 


c  \ 

eKM+P 

1-2-— I 

-hQ-X)" 

p_ 

l  P  J 

\  Po  J 

A-1 


+  t5(  1-2) 


+  A 


f  1  2  ^ 

r  r,  "i 

(  1  2  ^ 

f  R2  ) 

- — 

exp 

+  B 

- “ 

exp 

{cop  R^j 

l  xp) 

{cop  R2J 

{  XP) 

+  (Q  +  e0)A  (4.1.10) 


’ii  . . 


' 

(  _  ) 

N  1 

(1-2)w+1 

— 

-0-2) 

> 

{Po) 

, 

Although  (4.1.10)  is  complicated,  it  is  in  a  convenient  form  for  differentiation  through  the 
use  of  the  quotient  rule.  We  also  note  that  (4. 1.10)  consists  of  a  sum  of  eight  terms,  i.e., 

P  =  (4.1.11) 

L*  ;=i 


so  we  may  use  linearity  and  differentiate  each  term  individually.  If  we  designate  a  non¬ 
conservative  variable  of  differentiation  as  q  ,  q  e  {p,  2 ,  e} ,  then  we  have  that 


dP 

dq 


1  8 


dD 


D  — —  -  qt - 

dq  dq 


(4.1.12) 
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Equation  (4.1.12)  presents  a  very  convenient  method  for  evaluating  pressure  derivatives. 
Below,  we  list  explicit  equations  required  in  evaluating  (4.1.12). 


<h-i;  ^=0;  §^  =  0; 

0/7  52  de 


P 


=  n.  5]h=PcL.  5%=_1.  ^2  _ Q 


/7p  —  1/^/  9Co  —  A-'?  —  o,  y 

/70  0/7  p  OX  de 
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0  X 
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(4.1.13) 

(4.1.14) 

(4.1.15) 

(4.1.16) 

(4.1.17) 

(4.1.18) 

(4.1.19) 
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7s  =(  1-*)W+1 

dll 


(  ^ 

P_ 

vA)  y 


M-1;  c8=%;  ^  =  ^(1-1)- 


^  3p  p0 


/  A 

P_ 

V  A)  y 


AM 


(4.1.20) 


dA 


=  1-(N  +  1) 


P(  1-A) 

A>  J 


5t7s 


=  0 


We  also  have  that 


3D_ _ 1_  SD__1. 

3/?  /?2(y  ’3/1  g  ’ 


(4.1.21) 


Clearly,  we  may  use  (4.1.12)  through  (4.1.21)  to  evaluate  the  pressure  derivatives 
required  by  the  eigen-space  decomposition  discussed  in  Section  3. 

4.2  Finite  Volume  Discretization 


Ultimately,  we  must  discretize  the  governing  equations  (2.1.1)  in  order  to 
numerically  solve  the  detonation  problem.  We  may  illustrate  the  discretization  procedure 
by  considering  a  simplified  fonn  of  (2.1.1),  i.e., 


3U  3F 

dt  dx 


(4.2.1) 


where  S  is  a  vector  containing  all  of  the  source  terms.  To  enact  the  finite  volume 
discretization,  we  integrate  (4.2.1)  in  1-D  space  as  follows 


(4.2.2) 


Moreover,  we  obtain 


(4.2.3) 


Since  the  limits  are  fixed  in  the  first  tenn  of  (4.2.3)  and  since  we  assume  that  U  is 
continuous  on  the  interval  {xi_V2,xM/2) ,  we  may  interchange  the  order  of  integration  and 
differentiation  to  find  that 

p\  *i+l  /  2  *f+l  /  2 

—  {vdx  +  ¥\XMn=  [S  dx  (4.2.4) 

Qf  J  *i'—l  /  2  J 

*1-1/2  *i-l/2 

By  observing  that  the  integral  in  the  first  term  is  taken  over  space,  we  may  evaluate  it  as 
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(4.2.5) 


■*7+1/2 

\vdx  =  \Ji{xi+V2-xi_V2) 

xi- 1/2 

where  Uf  is  the  average  of  U  =  U(x,/)  taken  over  space  in  the  interval  \xMI2,xi_^l2\. 
This  interval  defines  cell  i  in  the  finite  volume  grid.  Because  of  the  integration,  observe 
that  IT  =  U,  (x)  .  If  we  also  apply  this  idea  to  the  source  term,  (4.2.4)  becomes 


dV 


dt 


(-*7+1/2  X/-1/2)  +  -f  r,  1/2  —  (X/+1/2  A-l/2) 


the  so-called  semi-discrete  form.  Hence, 


ClU; 


L  +■ 


dt  Xi+V2  Xi- 1/2 


(^i+1/2  f'M/2)  -  S,- 


(4.2.6) 


(4.2.7) 


The  values  of  F  used  in  (4.2.7)  are  evaluated  at  cell  interfaces  (natural  locations  for 
possible  discontinuities  in  Euler  solutions).  As  a  result,  at  each  interface,  F  is  evaluated 
as  a  numerical  flux  through  the  use  of  an  upwind  discretization  scheme  based  on  the 

values  of  U,  defined  at  the  cell  centers.  The  upwind  scheme,  described  later  in 
Subsection  4.4,  makes  use  of  the  theory  developed  in  Section  3. 


4.3  Temporal  Discretization 

The  semi-discrete  form  (4.2.7)  offers  certain  numerical  advantages  (or 

disadvantages,  depending  on  your  point  of  view).  This  form  effectively  decouples  the 

temporal  discretization  scheme  from  the  spatial  discretization.  As  a  result,  we  are  free  to 

choose  different  methods  for  each  discretization.  On  the  other  hand,  one  may  argue  that  it 

is  unwise  to  decouple  the  time  and  space  schemes.  Why?  Our  shock-capturing  scheme 

12 

fundamentally  relies  on  solutions  of  the  Riemann  problem  and  on  characteristics. 
Characteristics  adjoin  the  time  and  space  coordinates  in  an  inextricable  manner,  so  in  the 
strictest  sense,  these  coordinates  cannot  be  decoupled.  This  effect  has  led  to  the  creation 
of  a  large  family  of  schemes  based  upon  Godunov’s  method  that  couple  the  time  and 
space  discretization.13  Although  we  do  not  disagree  with  these  ideas,  our  development  is 
evolutionary,  so  it  is  very  important  that  we  understand  our  space  scheme  at  a 
fundamental  level.  For  these  reasons,  we  will  use  the  decoupled  approach  involving  what 
is  perhaps  the  simplest,  explicit  temporal  discretization  method.  Let  us  recall  (4.2.7)  and 
discretize  the  time  derivative  with  a  simple  forward  difference.  The  current  time  level  is 
indicated  by  the  superscript  n  . 


up-u" 

At 


+  -f(Kvz-Kvz)  =  s; 
Ar, 


(4.3.1) 


where  At  =  t"']  -t"  is  the  numerical  time-step,  and  Ax;  =  xMn  —  xi_x/2  is  the  spatial 
stepsize.  Note  that  (4.3.1)  represents  a  fully  explicit  method;  by  rearranging,  we  obtain 
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(4.3.2) 


U"+1  =U"  +A  t 


sr- 


/+1/2 


-F,-i 


n 


Ay 


Basically,  equation  (4.3.2)  implements  the  Euler  time  integration  method.14  The  only 
numerical  stability  control  we  place  on  (4.3.2)  involves  a  restriction  on  the  time-step  At . 
This  restriction  is  enforced  through  a  Courant-Friedrichs-Lewy  (CFL)  criterion.  We 
apply  a  factor  of  0.5  to  the  new  predicted  time-step  given  by 


A  tpred  = 


min 

\<i<i  max 


Ay, 


\u.\  +  \a; 


(4.3.3) 


4.4  The  Numerical  Flux 

As  we  mentioned  earlier,  the  flux  vector  F  defined  at  each  interface  must  be 
evaluated  via  an  upwind  method  in  order  to  facilitate  the  automatic  capturing  of  shock 
waves  without  numerical  oscillations.  Our  upwind  method  of  choice  is  Roe’s  flux 
difference  splitting  scheme.  ~  To  promote  notational  clarity,  let  us  designate  the 
numerical  flux  vector  by  the  symbol  f  while  retaining  the  symbol  F  for  the  regular  flux 
vector  (2.1.3)  defined  by  the  reactive  Euler  equations.  Roe’s  numerical  flux  vector  is 
simply  stated  below.11 


f  =  ^(T1+r'«-pll<u*-ut» 


(4.3.4) 


where  A  is  the  flux  Jacobian  matrix  defined  by  (3.3.23)  and  evaluated  at  the  interface  in 


R 


i+1 


i+1/2 

Figure  1.  Interface  Notation 

question.  The  (~)  notation  indicates  that  this  evaluation  is  conducted  with  the  use  of  Roe- 
averaged  variables.  The  designations  F  and  R  are  best  explained  by  referring  to  Figure  1 . 
The  subscript  F  or  R  designates  that  the  quantity  is  defined  just  to  left  or  right  of  the 
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interface,  respectively.  In  Figure  1,  the  interface  is  located  at  xi+U2  between  cell  i  and 

cell  i  + 1 .  Why  would  the  left  and  right  interface  values  of  some  property  differ?  The 
answer  is  very  simple.  Remember  that  we  stated  earlier  that  our  method  involves 
solutions  of  the  Riemann  problem.  These  solutions  admit  discontinuities,  e.g.,  shock 
waves.  Hence,  by  the  nature  of  a  discontinuity,  the  properties  taken  to  the  left  and  the 
right  of  an  interface  differ.  In  the  simplest  view,  we  can  say  that  the  properties  to  the  left 
of  the  interface  taken  on  the  values  defined  in  cell  i ;  it  follows  that  the  properties  to  the 
right  of  the  interface  take  on  the  values  defined  in  cell  i  + 1 .  This  means  of  selecting  the 
left  and  right  interface  values  renders  first-order  accuracy  on  unifonn  meshes.  There  are 
other  ways  to  define  these  upwind  values.  A  higher  order  method  is  discussed  in  a  later 
subsection.  Our  Roe  averages  are  computed  from  these  upwind  (L  and  R)  variables. 

The  Roe  average  constitutes  the  physically  correct  representation  of  an  average  at 
a  discontinuity  conforming  to  the  basic  ideas  of  flux  difference  splitting. 15  A 
mathematically  lengthy  derivation  is  required  to  produce  Roe’s  fonnulas,  so  we  merely 
state  the  results.10 


P  ~  yl  Pl  Pr 

(4.3.5) 

~  ul  yl  Pl  ur  yl  Pr 

VA+VA 

(4.3.6) 

H  H  l  yl Pl  +  Hr  yl  Pr 
yl  Pl  +  yl  Pr 

(4.3.7) 

~  eL  yl  Pl  +  eR  yl  Pr 
€  = 

yl  Pl  +  a/  Pr 

(4.3.8) 

l  _^l  yl  Pl  +  yl  Pr 

yl  Pl  +  yl  Pr 

(4.3.9) 

(„  \  \ 

P  =  p  H  —  e — u 2 

l  2  J 

(4.3.10) 

~2  n  PPe 

a-=Pp+  ' 

(4.3.11) 

One  may  note  that  (3.3.20)  through  (3.3.22),  (3.3.25)  and  (4.3.11)  require  Roe-averaged 
pressure  derivatives.  Recall  that  explicit  formulas  for  these  derivatives  are  presented  in 
(4.1.12)  through  (4.1.20).  The  derivatives  are  presented  in  tenns  of  the  primitive 
variables,  so  we  claim  that  Roe-averaged  values  of  the  pressure  derivatives  may  be 
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obtained  by  simply  evaluating  these  fonnulas  for  the  Roe-averaged  variables  presented  in 
(4.3.5)  through  (4.3. 10).  In  practice,  this  procedure  seems  to  work  well. 

We  may  now  address  the  practical  evaluation  of  the  numerical  flux  vector  as  it  is 
defined  in  (4.3.4).  The  vectors  FL  and  FA,  are  the  standard  Euler  flux  vectors  (2.1.3) 
evaluated  for  the  upwind  conservative  variables  UL  and  UA,  (or  primitive  variables  qA, 
and  qL),  respectively.  The  remaining  tenn 

a|(U*-UJ  (4.3.12) 

is  denoted  as  the  numerical  viscosity  expression.  The  difference  between  the  conservative 
variables  left  and  right  of  the  interface  may  be  easily  evaluated  through  the  use  of  (2.1.2). 

A  may  be  evaluated  as  follows. 

A  =RAL  (4.3.13) 

where  the  (~)  notation  indicates  that  all  of  the  entries  in  the  matrices  are  calculated  with 
the  use  of  averaged  variables.  The  matrix  A  is  created  by  taking  the  absolute  value  of 

each  element  of  A,  the  diagonal  matrix  of  eigenvalues.  Finally,  (4.3.12)  is  computed  by 
a  series  of  simple  matrix -matrix  and  matrix -vector  multiplications;  (4.3.4)  is  easily 
evaluated  by  using  vectors  sums. 

4.5  A  Higher-Order  Scheme 

The  scheme  described  in  the  preceding  subsection  is  only  accurate  to  the  first 
order,  and  it  is  highly  dissipative,  a  detriment  to  the  sharp  resolution  of  detonation  waves. 
In  this  subsection,  we  briefly  describe  an  enhancement  to  the  first  order  scheme  that  is 
third-order  accurate  on  uniform  grids.  As  you  may  have  concluded,  the  left  and  right 
interface  values  are  constructed  from  the  cell-center  values  to  the  left  and  right  of  the 
interface,  respectively.  To  increase  the  order  of  accuracy  for  the  scheme,  we  instead 
reconstruct  the  interface  values  using  interpolating  polynomials  involving  more  than  one 
cell-center  value.  One  way  to  apply  this  idea  is  through  the  use  of  a  Monotone  Upwind 
Scheme  for  Conservation  Laws  (MUSCL).  The  equations  for  the  left  and  right  interface 
variables  are  provided  below  for  the  interface  located  at  i— 1/2.  Consider  the  primitive 
variable  cj  ,  q  e  {p,  u,  P,  A}. 


<c,  =  Wi + t  (1-^)0(fi)(wi-^-2)+(1+^)°  —  k  -wi) 


(4.4.1) 


where  k-  1/3  to  achieve  third-order  accuracy,  and 
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Qi-Qi- 1 


(4.4.2) 


7  ;-l  C1  i-2 

O  is  a  function  designed  to  serve  as  a  non-limiter  limiter.  In  every  case,  our  interpolated 
data  must  be  monotone;  otherwise,  the  interpolation  procedure  will  result  in  the 
formation  of  non-physical  oscillations  in  the  numerical  solution.12  The  nonlinear  limiter 
is  designed  to  maintain  the  monotonicity  of  smooth  sections  of  data  when  interpolated  to 
high  order.  We  have  chosen  the  Van  Albada  limiter  for  use  in  this  problem,  i.e., 

^  (4-4.3) 

1  +  r 


The  right  interface  variable  is  given  by 


Qr  —  V 


(1  -  k)  ® (rR )(qM  -  q, )  +  (1  +  k)  ® 


\rRj 


1) 


(4.4.4) 


For  this  expression,  the  ratio  used  by  the  limiter  is  defined  as 


_  V  i 
Qi+i-Qi 


(4.4.5) 


Equations  (4.4.1)  through  (4.4.5)  cannot  be  implemented  without  due  cognizance.  The 
left  interpolant  involves  cell-center  values  located  at  i-2,  i- 1  and  i .  As  a  result,  we 
must  ensure  that 


(V  -  V  i )  (V  i  -  q>- 2 )  >  0  (4.4.6) 

Otherwise,  the  cell-center  data  is  non-monotone,  and  the  interface  values  must  be  set  to 
the  first-order  values 

9l  ~  (4.4.7) 

c1r  =  cli 

in  order  to  properly  smooth  the  solution.  For  the  right  interpolant,  we  must  ensure  that 

(q,  ~  q>- 1 )  (qi+ 1  ~qt)>  0  (4.4.8) 

or  we  must  use  the  first-order  interpolation  values  (4.4.7).  In  addition,  after  the  criteria 
(4.4.6)  and  (4.4.8)  are  satisfied,  we  are  required  to  limit  on  the  ratios  (4.4.2)  and  (4.4.5). 
Based  on  the  data,  these  ratios  may  become  undefined,  so  the  limiter  function  (4.4.3) 
must  be  modified  ensure  that  its  value  always  remains  finite.  If  this  interpolation  strategy 
is  used  properly,  the  Roe  algorithm  becomes  a  high-resolution  flux  difference  splitting 
scheme. 
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4.6  Boundary  Conditions 

In  most  cases,  we  cannot  solve  partial  differential  equations  without  applying 
boundary  conditions.  Even  for  our  simple  detonation  problem  cast  in  one  dimension,  we 
must  apply  boundary  conditions  at  x  =  0  (the  center  of  the  sphere)  and  at  x  =  xMAX  (the 
outer  surface  of  the  sphere).  At  the  center  of  the  sphere,  we  enforce  fully  reflective 
boundary  conditions  through  the  use  of  a  ghost  cell  installed  at  i  =  0 ,  i.e., 

Po  =  P\ 

VI  Q  VI  -j 

P0=  P ,  (4.5.1) 

K  =  4 

eo  =  e\ 

We  have  assumed  that  the  first  flow  field  cell  adjacent  to  this  boundary  has  the  index 
i  =  l. 

At  the  outer  surface  of  the  sphere,  we  apply  extrapolated  boundary  conditions  to 
mimic  a  supersonic  outflow.  We  implement  this  condition  by  installing  a  ghost  cell  at 
i  =  /max  .  We  set  conditions  in  this  cell  as  follows. 

PlMAX  —  PlMAX-1 
WIMAX  —  WIMAX-1 

^TMAX  =  ^IMAX-l  (4.5.2) 

^IMAX  —  ^IMAX-1 
^IMAX  —  ^IMAX-1 

Boundary  conditions  (4.5.1)  and  (4.5.2)  function  well  for  the  detonation  of  a  finite 
spherical  mass  of  HMX. 
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5  PARTICLE  MOTION 


In  this  section,  we  extend  our  discussion  beyond  the  application  of  numerical 
detonation  literature  cited  thus  far.  Given  the  level  of  interest  in  Multiphase  Blast 
Explosives  (MBX),  it  is  desirable  to  incorporate  solid  particles  into  our  detonation 
programming.  This  effort  is  new,  so  our  treatment  of  solid  particles  is  limited,  to  a  certain 
extent.  Still,  our  particles  have  realistic  mass  and  finite  radii.  They  are  driven  by  the 
detonation  through  the  use  of  Lagrangian  laws  of  motion.  Our  particle  algorithms  have 
only  three  major  limitations: 

(i)  The  particle  collection  exists  in  the  diffuse  limit.  Particles  are  assumed  not  to 
interact  with  one  another. 

(ii)  Particles  are  assumed  to  exist  as  rigid  spheres.  The  do  not  deform  or  change 
phase  during  the  detonation  event. 

(iii)  This  model  is  restricted  to  one  dimension.  We  can  only  establish  initial 
particle  positions  along  a  single  ray. 

Based  on  these  assumptions,  we  can  investigate  the  efficacy  of  this  model  in  predicting 
the  post-detonation  conditions  for  a  mass  of  solid  HMX  loaded  with  particles. 

5,1  Coupling  Terms 

We  may  now  discuss  the  coupling  terms  (source  terms)  for  particles  presented  in 
equations  (2.1.1)  and  (2.1.6).  Fs  and  Qs  have  relatively  simple  descriptions.  Fs 
represents  the  transfer  of  momentum  between  the  gas  phase  and  the  particle  phase  while 
Q s  represents  the  similar  transfer  of  thermal  energy.  For  spherical  particles,  these  terms 
may  be  written  in  a  simple  form.6  Assume  that  the  total  number  of  particles  is  N  . 


77  4  3  dup 

Fs=-L^Pprp 

p= i 


3  '  p  p  dt 


p= i 


(5.1.1) 

(5.1.2) 


where  p  ,  rp  and  up  are  the  solid  density,  radius  and  velocity  of  the  pth  particle, 
respectively.  Therefore,  dup/dt  is  the  acceleration  of  the  pth  particle.  Also,  T  is  the 
temperature  of  the  gas  phase  at  the  surface  of  the  particle,  and  T p  is  the  particle 

temperature.  Actually,  T  is  the  Favre-filtered  temperature;  this  filtering  operation  is  used 
to  take  the  presence  of  turbulence  into  account.  Our  simulation  is  non-viscous,  so  we 
simply  set  T  equal  to  the  gas  phase  temperature  T  .  The  parameter  hp  is  the  heat  transfer 
coefficient  that  governs  the  transfer  of  thennal  energy  at  the  particle/fluid  interface.  In 
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general,  hp  is  experimentally  detennined.  By  specifying  (5.1.1)  and  (5.1.2),  we  can 

accurately  describe  the  coupling  between  the  gas  and  particulate  phases.  Of  course,  these 
equations  only  apply  to  particles  of  fixed  mass.  Additional  tenns  (including  mass 
conservation)  must  be  specified  for  particles  that  react  with  the  gas  phase. 


5.2  Particle  Laws  of  Motion 

The  detonation  physics  algorithms  incorporate  discrete,  finite-mass  particles,  so 
we  apply  Lagrangian  equations  for  tracking  the  movement  of  particles.  Let  xp  designate 

the  radial  coordinate  of  the  ptb  particle.  Then  we  have  that 


dxp 

clt 


=  up 


(5.2.1) 


The  particle  velocity  up  must  be  detennined  from  the  evolution  equation  given  by  a 

model.  We  have  two  alternatives  for  this  model;  the  first  is  called  the  “Spray  Model” 
which  may  be  described  as  follows.6 


chip 

dt 


3  CDfi  Rcp 
1 6  Pp  rp 


( u-up ) 


(5.2.2) 


where  the  particle  Reynolds  number  Re  is  defined  as 


Re 


P 


2r„P  | 

- \U~Un 

P  ' 


(5.2.3) 


The  drag  coefficient  for  the  particle  CD  is  conveyed  by  the  “Spray  Drag  Law”,  i.e., 


CD=\ 


24 


Re, 


2/3  y 


0.44 


Rep  <1000 
Rep  >1000 


(5.2.4) 


p ,  //  and  u  are  the  density,  dynamic  viscosity  and  velocity  of  the  gas  phase  in  the 
vicinity  of  the  particle.  This  model  is  not  appropriate  for  detonation  problems,  but  it  still 
serves  well  for  testing.  For  the  problem  of  a  detonation  with  solid  inclusions,  we  apply  a 
high  speed  gas  flow  model  originally  developed  for  solid  rocket  motors. 
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The  high  speed  gas  flow  model  was  developed  for  the  multiphase  flow  field 
created  by  the  burn  of  porous,  powdered  explosive  material.16  In  this  case,  the  particle 
acceleration  is  given  by 


chip 

dt 


n  d2pCDp 
8  mp 


U-Up  ( u 


(5.2.5) 


In  order  to  maintain  our  notation  consistent  with  the  literature,  (5.2.5)  is  written  in  terms 
of  the  particle  diameter  d  instead  of  the  radius.  Also,  mp  is  the  mass  of  the  p'h  particle. 
This  high  speed  drag  law  provides  the  drag  coefficient  through  a  more  complicated 
calculation.  First,  we  calculate  a  “Mach-zero”  drag  coefficient,  CD0,  i.e., 


Ci 


(0.45-a2)C1  +  (a2  -0.08 )C2 


0.37 


a2  <  0.08 


0.08  <a2<  0.45 


C0 


a2  >  0.45 


(5.2.6) 


where  Re  is  calculated  by  using  (5.2.3),  and 


C,= 

c2  = 


24  4.4 


+  0.42 


3  a. 


1.75  + 


150a2 
a\  RS  , 


(5.2.7) 

(5.2.8) 


In  (5.2.6)  and  (5.22.8),  we  have  introduced  two  new  parameters  a]  and  a2 ;  they  are  the 
volume  concentrations  of  the  gas  and  particle  phases,  respectively.  These  parameters 
require  interpretation  when  considering  the  detonation  problem.  At  the  outset  of  the 
problem,  the  solid  explosive  has  not  been  detonated,  so  there  is  no  gas  phase  at  this  point. 
The  best  course  of  action  is  to  compute  the  initial  values  of  ax  and  a2  based  upon  the 
volume  of  the  solid  explosive  and  the  volume  of  particles.  Since  we  are  not  simulating 
details  of  the  shock  interaction  with  metal  particles,  we  calculate  al  and  a2  on  this  basis 
of  the  initial  calculation  and  maintain  them  fixed  for  the  duration  of  the  detonation.  We 

17. 

must  then  calculate  a  final  value  of  CD  based  on  a  Mach  correction.  This  correction 
exists  due  to  the  natural  variation  in  the  drag  coefficient  with  Mach  number.  If  we  do  not 
wish  to  implement  a  drag  correction,  then  we  set  CD  =  CD0 ;  otherwise  the  corrected  value 
of  CD  may  be  calculated  from 
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(  (  o  427  Yl 

CD=CD0  1  +  exp  ,  (5.2.9) 

V  v  M  JJ 

where 

\u  —un\ 

M  =  l- - (5.2.10) 

a 

By  using  the  particle  velocities  provided  by  (5.2.2)  though  (5.2.4)  or  (5.2.5)  through 
(5.2.10),  we  may  integrate  (5.2.1)  to  determine  the  track  of  each  particle  through  space 
during  the  detonation. 
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6  RESULTS 


From  the  start  of  this  effort,  several  versions  of  our  current  numerical  detonation 
computer  code  have  been  developed  by  the  author.  The  purpose  of  this  section  is  to 
present  some  of  the  results  produced  for  typical  problems.  Specifically,  we  discuss  three 
results.  The  first  set  of  results  is  intended  to  show  that  our  detonation  program  is 
functioning  properly  and  producing  physically  correct  solutions.  In  a  second  calculation, 
we  address  the  numerical  detonation  of  a  spherical  mass  of  pure  HMX.  For  this  problem, 
we  have  computed  results  by  using  both  the  Hayes-I  and  Hayes-II  equations  of  state  for 
the  solid  explosive  combined  with  the  JWL  EOS  for  the  detonation  products.  Finally,  we 
discuss  the  results  for  the  detonation  of  a  spherical  mass  of  HMX  loaded  with  steel 
particles. 

6.1  Simple  Plane  Wave  Detonation 

This  test  problem,  described  in  Reference  2,  is  used  to  show  whether  or  not  the 
flux  difference  splitting  scheme  is  working  properly.  In  this  case,  we  endeavor  to  solve  a 
Deflagration  to  Detonation  Transition  (DDT)  problem  in  one  dimension.  Both  the 
explosive  and  the  detonation  products  are  modeled  by  using  the  calorically  perfect  gas 
EOS.  The  associated  mixture  EOS  is  given  as 

e  =  —^—-QX  (6.1.1) 

p(r- !) 


As  discussed  in  Section  4,  we  apply  fully  reflective  boundary  conditions  at  x  =  0  and 
extrapolation  conditions  at  x  =  xMAX.  For  this  problem,  we  use  the  reaction  rate 
expression 


r  =  k(  1-T)exp 


PJ 


(6.1.2) 


where  (6.1.2)  is  in  Arrhenius  fonn;  k  is  the  reaction  rate  constant,  and  Ea  is  a  parameter 
that  behaves  like  an  activation  energy.  The  one-dimensional  domain  is  defined  in 
0<x<12.  Also,  we  have  that  £’u=10;  <2  =  50;  y  =  1.4,  and  k  =  7.  The  problem  is 
initialized  with  u  =  0  ;  P  =  0 ,  and  A  =  0  everywhere.  The  initial  density  distribution  is 
given  by 

/?(*)=  Q  V  2,,  0<x<12.  (6.1.3) 

1  +  3exp(-x  ) 


This  density  distribution  initiates  the  reaction  in  the  region  near  x  =  0  by  boosting  the 
reaction  rate  term. 
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Density 


Figure  2.  Problem  1  Detonation  Field  Density,  Time  =  3.0 


Vel  ocity 


Figure  3.  Problem  1  Detonation  Field  Velocity,  Time  =  3.0 
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Figure  4.  Problem  1  Detonation  Field  Pressure,  Time  =  3.0 
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Figure  5.  Problem  1  Detonation  Field  Reaction  Progress  Variable,  Time  =  3.0 

This  problem  does  not  possess  an  “exact”  solution,  but  Xu  et  al.  have  obtained  a 
fully  converged  numerical  solution  using  a  mesh  consisting  on  3200  cells.2  This  problem 
provides  an  excellent  test  detonation  physics  algorithms.  Accordingly,  we  have  generated 
three  numerical  solutions  on  grids  comprised  of  200,  800  and  3200  cells,  respectively. 
The  numerical  solutions  for  density,  velocity,  pressure  and  the  reaction  progress  variable 
are  provided  in  Figures  2  through  5,  respectively,  at  the  dimensionless  time  3.0.  In  each 
figure,  solution  plots  are  color-coded  to  correspond  to  the  mesh  used.  The  behavior 
shown  in  each  plot  agrees  quite  well  with  archived  plots.  We  have  observed  only  one 
anomaly  in  our  solutions.  Strangely  enough,  on  the  mesh  consisting  of  only  200  cells, 
there  are  noticeable  oscillations  in  the  reaction  progress  variable.  These  oscillations 
dissipate  with  increasing  mesh  density.  The  explanation  for  this  behavior  is  not 
immediately  evident.  In  some  of  our  solutions,  the  reaction  progress  variable  has  been 
observed  to  hunt  between  the  solid  and  gaseous  equations  of  state.  In  fact,  this  variable  is 
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Reaction  Progress 


very  sensitive  and  couples  strongly  to  the  reaction  rate.  We  apply  no  post-solution 
filtering  to  this  variable.  Secondly,  we  are  using  a  weak  time  integration  scheme  with 
poor  numerical  stability  perfonnance.  The  oscillations  become  less  prevalent  with 
increasing  grid  density,  so  the  space  scheme  may  be  compensating  for  the  time  scheme. 
This  phenomenon  bears  further  investigation  as  this  work  continues.  We  will  also  re¬ 
examine  the  nonlinear  limiter  coding.  Nevertheless,  our  converged  solution  agrees  well 
with  the  converged  archival  solution." 

6.2  Detonation  of  Pure  HMX 

This  problem  is  intended  to  demonstrate  our  computer  code’s  capability  for 
simulating  the  detonation  of  a  sphere  of  pure  HMX.  This  problem  permits  a  test  of  our 
discretization  of  the  geometric  source  term  found  in  the  reactive  Euler  equations  (2.1.1) 
and  (2.1.4).  It  also  represents  our  first  attempt  at  capturing  the  physics  of  a  realistic 
detonation  event.  In  this  case,  we  address  the  detonation  of  sphere  of  solid  HMX  with  a 
radius  of  4.5  cm.  The  radius  of  the  sphere  is  divided  into  800  cells.  Figure  6  shows  the 
density,  velocity,  pressure  and  reaction  progress  variables  for  the  numerical  solution  at 
three  microseconds  (ps)  detonation  elapsed  time.  As  you  can  see,  the  Von  Neumann 
spike  is  clearly  resolved  in  this  solution  as  is  the  Taylor  wave.  Moreover,  the  Chapman- 
Jouquet  pressure  is  captured  at  the  experimentally  obtained  value  of  42  GPa.  Also,  the 
numerical  detonation  velocity  has  a  value  of  1.02  cm/ps  which  is  very  close  to  the 
experimentally  obtained  value  of  0.911  cm/ps.21  Of  course,  the  experimental  value  is 
generally  taken  from  tests  that  mimic  plane  wave  detonation  conditions.  As  a  result,  we 
expect  to  calculate  a  different  value  for  the  spherical  detonation  problem.  Overall,  the 
results  agree  very  closely  with  the  archival  data.  We  have  also  solved  this  same  problem 
by  using  the  Hayes-II/JWL  mixture  EOS.  The  results  of  this  analysis  are  given  in  Figure 
7.  It  is  interesting  to  observe  that  the  Taylor  wave  is  captured  in  this  solution  even  more 
smoothly  than  it  was  in  the  preceding  case.  The  more  complex  Hayes-II  EOS  may 
actually  offer  greater  stability  when  used  in  the  mixture  EOS.  This  numerical  solution 
also  offers  excellent  comparisons  with  the  Chapman-Jouquet  pressure  and  detonation 
velocity  for  HMX.  Both  mixture  equations  of  state  show  that  the  detonation  reaction 
occurs  in  a  nearly  instantaneous  manner.  As  you  can  see,  the  reaction  progress  variable 
changes  in  a  nearly  discontinuous  manner  at  the  detonation  front.  In  either  case,  our 
computer  programming  captures  the  appropriate  physics  for  the  detonation,  and  it  renders 
a  wide  array  of  physical  data  (far  more  than  is  shown  here). 
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0«w«  tftgmOl  Woo w  i  vacs  I 


Figure  6.  Numerical  detonation  solution  Hayes-I/JWL  in  HMX  at  3  ns.  Horizontal  axis  is  distance  in 
meters. 
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Figure  7.  Numerical  detonation  solution  Hayes-I/JWL  in  HMX  at  3  fis.  Horizontal  axis  is  distance  in 
meters. 
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6.3  Detonation  of  HMX  Containing  Metal  Particles 

This  test  case  is  the  final  detonation  problem  addressed  by  this  report.  We 
consider  the  detonation  of  a  spherical  mass  of  HMX  loaded  with  a  radial  distribution  of 
steel  particles.  The  mass  of  the  HMX  sphere  remains  the  same  as  is  used  for  the 
preceding  problem,  and  we  still  have  800  finite  volume  cells  defined  along  the  charge 
radius.  For  this  example,  we  have  placed  ten  particles,  at  uniform  spacing,  along  the 
charge  radius.  The  particles  each  have  a  radius  of  463  pm  and  a  material  density  of  7860 
kg/m  .  We  assume  the  gas  viscosity  has  a  value  of  1.7x10'  kg/(m.s).  Furthermore,  in  this 
simulation  study,  we  have  applied  the  high  speed  flow  drag  law.  The  results  for  particle 
locations  are  presented  in  Figure  8  while  the  plot  of  particle  velocities  is  given  in  Figure 
9.  The  particle  tracks  shown  in  Figure  8  clearly  indicate  the  passage  of  the  detonation 
wave.  For  particles  farther  away  from  the  charge  center,  the  particle  tracks  show  changes 
in  slope  at  progressively  larger  times.  The  sudden  change  in  track  slope  concurs  with  the 
nearly  discontinuous  change  seen  in  the  particle  velocity  traces  shown  in  Figure  9.  Also, 
in  Figure  9,  the  effect  of  the  drag  law  can  clearly  be  seen  as  the  particle  velocities  rise 
rapidly  in  the  wake  of  the  detonation  wave  then  fall  quickly  under  the  action  of  drag  in 
the  region  behind  the  wave.  We  have  also  applied  the  Mach  correction  to  the  rocket  drag 
law.  In  the  velocity  trace  for  the  particle  closest  to  the  charge  center,  we  can  see  the 
velocity  begin  to  level  off  at  4.5  ps.  Available  data  indicates  that  the  calculated  terminal 
velocity  at  or  near  375  m/s  is  an  acceptable  value.  This  simulation  does  not  include 
thermal  effects  since  we  are  still  in  the  process  of  completing  our  detonation  products 
EOS. 


Particle  Location  (m) 


Time  (s) 

Figure  8.  Radial  locations  for  steel  particles  embedded  in  a  mass  of  detonating  HMX 
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Particle  Velocity  (m/s) 
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7  CONCLUSIONS 


In  this  report,  we  have  presented  the  governing  equations  for  the  direct  numerical 
simulation  of  the  detonation  of  a  solid  explosive  material.  Proper  equations  of  state  have 
been  discussed  for  both  the  solid  explosive  material  and  for  the  gaseous  detonation 
products.  From  these  equations  of  state,  we  have  developed  a  mixture  equation  of  state 
relating  the  specific  internal  energy  for  the  detonation  to  the  thermodynamic  pressure. 
The  resulting  computer  program  has  been  tested  on  an  archival  detonation  problem  for 
the  purpose  of  comparison.  We  have  presented  results  for  the  detonation  of  a  spherical 
mass  of  pure  HMX. 

More  importantly,  we  have  incorporated  particle  tracking  algorithms  within  the 
programming.  As  a  result,  the  code  can  now  explosively  drive  particles  under  the  action 
of  a  detonation  wave  with  coupling  to  a  drag  law.  This  mechanism  allows  the  code  to 
simulate  the  detonation  of  a  Multiphase  Blast  Explosive  in  the  diffuse  limit  of  particle 
loading.  We  have  built  drag  laws  for  both  spray  and  high  speed  gas  flow  drag  law  into  the 
code.  For  a  test  problem,  we  have  simulated  the  detonation  of  a  mass  of  HMX  loaded 
with  a  radial  distribution  of  steel  particles.  The  trend  in  post-detonation  velocities  of  these 
particles  meet  our  expectations. 

8  RECOMMENDATIONS 

During  the  months  ahead,  detonation  physics  algorithms  are  scheduled  for 
implementation  in  LESLIE3D.  The  development  of  the  present  work  has  been  a  learning 
experience  accompanied  by  a  large  number  of  difficulties,  especially  in  the 
implementation  of  Roe’s  flux  difference  splitting  scheme.  A  first  recommendation  is  that 
the  HLL  family  of  schemes  be  used  instead.  These  schemes  are  more  robust  and  do  not 
require  the  use  of  pressure  derivatives.  Also,  these  schemes  already  operate  well  inside  of 
LESLIE3D.  The  detonation  physics  solver  will  also  benefit  from  the  interface  tracking 
scheme  already  coded  into  LESLIE3D.  Clearly,  the  governing  equation  differ  at  the 
interface  between  the  condensed  explosive  and  the  surrounding  gas  field.  This  situation 
necessitates  an  interface  to  maintain  code  stability. 

The  detonation  physics  algorithms  discussed  here  must  be  adapted  for  curvilinear 
coordinates  in  three  dimensions.  For  HLL  flux  forms,  this  process  should  not  be  difficult. 
The  author  has  already  done  some  work  in  this  area.  However,  the  pressure  and  specific 
volume  (or  density)  closures  associated  with  the  mixture  equation  of  state  do  require 
attention.  The  Gas-Interpolated  Stewart-Prasad-Asay  (GISPA)  method  requires  these 
closures  to  address  the  multiphase  physics  of  detonation.  There  is  no  unique  set  of 
closures  available  for  this  process,  but  the  chosen  closures  must  be  carefull  accomplished. 
Some  difficulty  has  been  encountered  in  the  use  of  the  specific  volume  closure  (due  to 
Xu),  and  this  difficulty  should  be  investigated  and  resolved. 

The  Hayes  equation  of  state  for  the  solid  explosive  is  an  older  relationship  that 
characterizes  very  few  explosives.  The  Mie-Gruneisen  equation  of  state  characterizes 
many  more  explosive  materials.  That  is  to  say,  there  is  data  available.  However,  the 
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mixture  equation  of  state  must  be  rederived  for  the  Mie-Gruneisen  formulation.  It  may  be 
combined  with  the  JWL  adiabat  for  the  detonation  products,  or  with  another  real  gas  state 
equation.  The  “Wide-Ranging”  detonation  equation  of  state  may  also  be  implemented.4 

Ultimately,  the  particle  phase  algorithms  discussed  here  must  be  rewritten  for 
dense  phase  fields.  The  detonation  of  a  condensed  explosive  with  solid  inclusions  is  a 
dense  phase  problem.  Also,  the  computer  program  is  currently  not  properly  written  even 
in  the  diffuse  limit  as  regards  the  nonhomogeneous  source  terms.  The  integration  scheme 
should  be  changed  to  reflect  the  use  of  Strang  splitting.1  That  is  to  say,  the  spatial 
integration  scheme  should  be  advanced  in  separate  step  from  the  nonhomogeneous  tenns. 
For  the  latter  step,  the  integration  should  be  conducted  in  the  temporal  manner  at  each 
grid  cell  just  like  an  initial  value  problem. 
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APPENDIX  A 
SOURCE  CODE 


Instructions: 

The  source  code  that  follows  has  been  developed  over  a  period  of  six  years,  but  in 
a  sporadic  manner,  as  time  has  pennitted.  FORTRAN  77  is  used  throughout  the  computer 
program,  and  an  in-line  coding  structure  has  been  used.  The  programming  is  designed  for 
research  and  is  thus  rather  crude.  The  initial  conditions  (shock-based  initiation)  are  all 
rigidly  coded.  Different  initiation  options  exist,  but  they  must  be  enabled  or  disabled  by 
commenting.  The  detonation  reaction  rate  laws  are  treated  in  the  same  way.  The  desired 
reaction  rate  law  must  be  commented  in  for  the  initial  conditions  and  for  the  first  and 
second  time  step  segments  of  the  solver.  The  calorically  perfect  gas  and  Jones -Wilkins- 
Lee  test  problems  are  also  activated  or  deactivated  by  commenting  in/out  code  segments. 

This  computer  program  is  written  for  standard  explosives  like  HMX  for  which  we 
have  plenty  of  data.  Especially  for  the  Hayes  equation  of  state,  a  great  deal  of  data  input 
is  required.  This  data  is  simply  entered  directly  into  the  source  code.  This  statement  is 
also  true  as  pertains  to  the  Jones-Wilkins-Lee  detonation  product  data  as  well  as  the 
particle  field  data.  This  code  functions  in  one  dimension  only:  Cartesian,  cylindrical  or 
spherical.  The  domain  boundaries  are  contained  between  xl  and  x2.  The  number  of  cells 
in  the  detonation  field  is  given  by  imax-1.  The  variable  NSTP  tells  the  code  how  many 
iterations  (time  steps)  to  execute  while  the  variable  NDMP  tells  the  code  how  many 
iterations  to  perform  between  dump  files.  The  variable  IRST  controls  code  execution. 
With  IRST  set  at  zero,  the  code  begins  with  the  coded  initial  conditions.  With  IRST  set  at 
one,  the  code  reads  the  restart.data  file  to  obtain  its  starting  conditions.  The  IEOS 
variable  switches  between  the  mixture  equations  of  state.  IEOS  equal  zero  sets  calorically 
perfect  gas  conditions.  IEOS  at  one  sets  JWL  conditions  while  IEOS  equal  2  or  3  sets  the 
Hayes-EJWL  and  Hayes-II/JWL  formulations.  The  reader  should  be  advised  that  the  pure 
JWL  option  does  not  work  well.  The  fault  of  this  equation  is  that  there  is  not  a  sufficient 
energy  separation  between  the  adiabats  to  result  in  detonation. 

This  detonation  physics  program  utilizes  a  number  of  flags  and  control  parameters 
in  order  to  stabilize  code  operation.  Some  of  these  parameters  set  tolerances  on  the 
variables  (like  the  reaction  progress  variable)  to  prevent  “hunting”.  Other  flags  control 
solution  progress.  For  instance,  internal  energy  updates  are  lagged  by  one  iteration  to 
keep  temperature  from  turning  negative.  It  is  also  important  to  observe  that  the  equations 
of  state  used  here  have  constant  specific  heat  formulations.  Over  time,  this  limitation 
should  be  lifted,  but  better  equation  of  state  data  is  required  to  do  so.  We  also  zero  the 
detonation  reaction  rate  in  the  far  field.  As  it  happens,  the  flux  scheme  will  erroneously 
allow  reaction  rate  to  creep  up  slowly  in  the  unreacted  explosive  mass.  This  effect  is 
damaging  to  the  solution  and  had  to  be  corrected. 

c*************  EZ1_MASTER  ************ 

(-I'k-k'k-k-k-k-k-k'k-k-k-k'k-k'k-k-k'k-k-k-k'k-k-k-k'k-k'k-k'k-k 

c  Program  for  1-D  detonation  test  problem 
c  Simple  coding  structure 
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c  Monotonicity  check  implemented  on  extrapolation 
c  Direct  adaptations  for  calorically  perfect  gas  and  JWL 

program  ezl  master 
implicit  none 

c  Parameter  statements 
integer  imax 

c  parameter  (imax  =  20001) 

parameter  (imax  =  2001) 

integer  npar 
parameter  (npar  =  1000) 

real*8  cl2 

parameter  (cl2  =  0.5d0) 
real*8  cl3 

parameter  (cl3  =  Id0/3d0) 
real*8  cl4 

parameter  (cl4  =  0.25d0) 
real*8  cl8 

parameter  (cl8  =  0.125d0) 
real*8  c23 

parameter  (c23  =  2d0/3d0) 
real*8  c43 

parameter  (c43  =  4d0/3d0) 
real*8  c316 

parameter  (c316  =  3d0/16d0) 
real*8  pi 

parameter  (pi  =  3 . 141592654d0) 

c  Variable  array  declarations 
c  File  I/O 

character*12  filex 
character*12  parex 

c  Debug  flags 

integer  idbgl 
integer  idbgf 
integer  idbgs 
integer  idbgp 

c  Control  flags 

integer  irst 
integer  ieos 
integer  igeo 
integer  irxn 
integer  ipar 
integer  idrg 
integer  imach 
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integer  iext 
integer  iav 
integer  ilim 
integer  imon 
integer  iefx 
integer  item 

c  Counters 

integer  i 
integer  n,nn,np 
integer  l,m 
integer  k 
integer  nstart 
integer  nstp 
integer  ndmp 
integer  nfil 

c  Gas  phase  data 
real* 8  pamb 
real*8  mu 

c  Calorically  perfect  EOS  data 
real*8  gamm 
real*8  garni 

c  JWL  EOS  data 
real*8  rO 
real*8  aj 
real*8  bj 
real*8  cj 
real*8  cjh 
real*8  rl 
real*8  r2 
real*8  wj 
real*8  pcj 

c  Hayes-I  EOS  data 
real*8  cvs 
real*8  gh 
real*8  hi 
real*8  nh 
real*8  rgas 
real*8  cvg 
real*8  cpg 
real*8  nhpl 
real*8  nhml 
real* 8  nhm2 
real*8  t3 
real*8  t4 
real*8  t5 
real*8  t7 
real*8  alfa 
real*8  beta 
real*8  thta 

c  Mixture  EOS  tolerances 
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real*8  ztoll 
real*8  ztol2 


c  Detonation  data 
real*8  qdetO 
real*8  eO 
real*8  eact 
real*8  rk 
real*8  rkl 
real*8  rk2 
real*8  pexp 
real*8  zexp 
real*8  thl 
real*8  th2 

real*8  rhl 
real*8  rh2 
real*8  rht 
real*8  rhti 
real*8  wrl 
real*8  wr2 
real*8  wrlr 
real*8  wr2r 

c  Grid/Timestep  control  data 
real*8  xl 
real*8  x2 
real*8  chx 
real*8  dx 
real*8  xc 
real*8  fct 
real*8  fctl 
real*8  fct2 

real*8  time 
real*8  tend 
real*8  dt 
real*8  dtO 
real*8  dtl 
real *8  dtmx 
real*8  cfl 
real*8  offs 

c  Derived  data 
real*8  et 
real*8  ra 
real*8  ra2 
real*8  ea 
real*8  za 
real*8  rz 
real*8  omz 
real*8  rxmin 

real*8  bot 
real*8  bot2 
real*8  botr 
real*8  botz 

Distribution  A.  Approved  for  public  release,  distribution  unlimited.  (96ABW-20 1 1  -0548) 

45 


real*  8 

dpdr 

real*  8 

dpde 

real*  8 

dpdz 

real*  8 

a2 

real*  8 

psgn 

real*  8 

kap 

real*  8 

eps 

real*  8 

epsm 

real*  8 

epsp 

real*  8 

off 

real*  8 

tmp 

real*  8 

rl ,  rr 

real*  8 

ul ,  ur 

real*  8 

pl,pr 

real*  8 

zl ,  zr 

real*  8 

el ,  er 

real*  8 

eel , eer 

real*  8 

hhl , hhr 

real*  8 

dqer, dqwr, dqir 

real*  8 

dqeu , dqwu , dqiu 

real*  8 

dqep, dqwp, dqip 

real*  8 

dqez, dqwz, dqiz 

real*  8 

denm 

real*  8 

dra, drb, drc, drd, dre 

real*  8 

dua, dub, due, dud, due 

real*  8 

dpa, dpb, dpc, dpd, dpe 

real*  8 

dza,dzb,dzc,dzd,dze 

real*  8 

rat 

real*  8 

phir 

real*  8 

phiu 

real*  8 

phip 

real*  8 

phiz 

real*  8 

phi 

real*  8 

vhi 

real*  8 

sqrl 

real*  8 

sqrr 

real*  8 

rsumi 

real*  8 

rav 

real*  8 

ri 

real*  8 

uav 

real*  8 

zav 

real*  8 

eav 

real*  8 

hav 

real*  8 

aav 

real*  8 

pav 

real*  8 

delr 

real*  8 

delv 

real*  8 

delp 

real*  8 

delz 
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real*8  detr 
real*8  pest 

c  Temperature  estimation  variables 


real*  8 

tkO 

real*  8 

dtkmx 

real*  8 

denmx 

real*  8 

numr 

real*  8 

eOcr 

real*  8 

eta 

real*  8 

rs 

real*  8 

rg 

real*  8 

del 

real*  8 

de2 

real*  8 

de3 

real*  8 

de4 

real*  8 

de5 

real*  8 

de6 

c  Particle  phase  data 


real*  8 

xpl 

real*  8 

xp2 

real*  8 

dxp 

real*  8 

rdp 

real*  8 

dip 

real*  8 

rop 

real*  8 

pep 

real*  8 

rep 

real*  8 

ppr 

real*  8 

tcon 

real*8 

crppr 

real*  8 

nup 

real*  8 

hp 

real*  8 

cdp 

real*  8 

pum 

real*  8 

pam 

real*  8 

delu 

real*  8 

adelu 

real*  8 

hevol 

real*  8 

pvol 

real*  8 

cvol 

real*  8 

pOmas 

real*  8 

pmass 

real*  8 

alf  1 

real*  8 

alf  2 

real*  8 

alf  2 1 

real*  8 

cdl 

real*  8 

cd2 

real*  8 

cdO 

real*  8 

mach 

real*  8 

dtp 

c  Array  declarations 
real* 8  x ( imax) 
real*8  r(0:imax) 
real*8  p(0:imax) 
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real*  8 

u ( 0 : imax) 

real*  8 

z ( 0 : imax) 

real*  8 

ei ( 0 : imax) 

real*  8 

a ( 0 : imax) 

real*  8 

rxr ( 0 : imax) 

real*  8 

c  (8) 

real*  8 

top ( 8 ) 

real*  8 

topr ( 8 ) 

real*  8 

topz ( 8 ) 

real*8 

rp ( 0 : imax) 

real*  8 

pp ( 0 : imax) 

real*  8 

up ( 0 : imax) 

real*  8 

zp ( 0 : imax) 

real*  8 

eip ( 0 : imax) 

real*8 

etp ( 0 : imax) 

real*  8 

ap ( 0 : imax) 

real*  8 

tk ( imax) 

real*  8 

dtk ( imax) 

real*  8 

zzl (imax) 

real*  8 

zzr (imax) 

real*  8 

qv ( imax, 4 ) 

real*  8 

qvp ( imax, 4 ) 

real*  8 

sg ( imax, 4 ) 

real*  8 

srx ( imax, 4 ) 

real*  8 

sp ( imax, 4 ) 

real*  8 

s ( imax, 4 ) 

real*  8 

aeg ( 4 ) 

real*  8 

evr (4,4) 

real*  8 

cwm ( 4 ) 

real*  8 

chkl (4,4) 

real*  8 

chk2 (4,4) 

real*  8 

dq  ( 4 ) 

real*  8 

vl  (4) 

real*  8 

vn  (4) 

real*  8 

fl  (4) 

real*  8 

fr  (4) 

real*  8 

f n ( imax, 4 ) 

real*8 

dqv ( 4 ) 

real*  8 

derv ( imax, 2 

c  Particle  arrays 

integer  peel (npar) 
real*8  px(npar) 
real*8  pu (npar) 
real*8  pa (npar) 
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real*8  pxp(npar) 
real*8  pup(npar) 
real*8  pq(npar) 
real*8  ptk(npar) 
real*8  ptkp(npar) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c 

Main 

Data  Entry  Section 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c  Grid  data 

xl 

=  OdO 

c 

x2 

=  200d0 

x2 

=  3 . 6d-2 

chx 

=  3 . 8d-2 

c  CPG 

EOS  data 

gamm 

=  1 . 4d0 

pamb 

=  101325d0 

rgas 

=  287d0 

c  Extrapolation  control  data 

kap 

=  Id0/3d0 

c 

kap  = 

-IdO 

eps 

=  ld-12 

c  EOS 

control  tolerances 

ztoll 

=  ld-2 

c 

ztoll 

=  OdO 

ztol2 

=  0 . 99d0 

c 

ztol2 

=  IdO 

c  HMX 

Hayes 

EOS  Data  (Xu) 

c 

rO 

=  1891d0 

c 

hi 

=  1 . 35dl0 

c 

CVS 

=  1 . 5d3 

c 

gh 

=  2 . Id3 

c 

nh 

=  9 . 8d0 

c 

tkO 

=  3d2 

c  HMX 

JWL  EOS  Data  (Zukas/Xu) 

c 

aj 

=  7 . 783dll 

c 

bj 

=  0 . 07071dll 

c 

cj 

=  0 . 00643dll 

c 

rl 

=  4 . 2d0 

c 

r2 

=  IdO 

c 

wj 

=  0 . 3d0 

c 

cvg 

=  (2 . 4d0  -  0 ,28d0*r0*ld-3  -  1.3d0)*ld3 

c  NM  : 

Hayes 

EOS  Data 

c 

rO 

=  1 . 13d3 

c 

hi 

=  1 . 32d9 

c 

CVS 

=  1 . 446d3 

c 

gh 

=  1 . 356d3 

Distribution  A.  Approved  for  public  release,  distribution  unlimited.  (96ABW-20 1 1  -0548) 

49 


c 

c 


nh 

tkO 


7 . 144d0 
2  93d0 


c 

NM 

jwl  : 

EOS 

Data 

c 

aj 

= 

209. 2d9 

c 

bj 

= 

5 . 689d9 

c 

cj 

= 

0 . 77d9 

c 

rl 

= 

4 . 4d0 

c 

r2 

= 

1 . 2d0 

c 

wj 

= 

0 . 3d0 

c 

cvg 

— 

1 . 3d3 

c 

RDX 

Hayes  EOS  Data 

rO 

= 

1 . 6d3 

hi 

= 

13d9 

CVS 

= 

1 . 163d3 

gh 

= 

1 . 356d3 

nh 

= 

6 . 3d0 

tkO 

= 

300d0 

c 

RDX 

JWL 

EOS 

1  Data 

aj 

= 

573 . 187d9 

bj 

= 

14 . 639d9 

cj 

= 

0 . 77d9 

rl 

= 

4 . 6d0 

r2 

= 

1 . 4d0 

wj 

= 

0 . 32d0 

cvg 

= 

1 .2d3 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Detonation  reaction  data 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c 

CPG  Test 

eact 

= 

lOdO 

rk 

= 

16 . 418d0 

c 

thl 

= 

OdO 

c 

th2 

= 

OdO 

c 

rxmin 

= 

rk*dexp (-eact) 

c 

qdetO 

25d0 

c 

HMX  Test 

c 

pcj 

- 

42d9 

c 

rkl 

= 

110d6 

c 

rk2 

= 

OdO 

c 

pexp 

= 

3 . 5d0 

c 

zexp 

= 

0 . 93d0 

c 

thl 

= 

OdO 

c 

th2 

= 

OdO 

c 

rxmin 

= 

rkl* ( (pamb/pcj ) **pexp) 

c 

qdetO 

= 

( 7 . 91d0  -  4 . 33d0* (r0*ld-3  -  1.3d0)**2 

c 

& 

-0 . 934d0* (r0*ld-3  -  1.3d0))*ld6 

c 

NM  Test 

c 

pcj 

= 

12 . 5d9 

c 

pexp 

= 

ldO 

c 

zexp 

= 

0 . 95d0 

c 

rkl 

= 

7 . 75dl0 
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c 

rk2 

= 

1 . 5dl2 

c 

thl 

= 

14500d0 

c 

th2 

= 

29700d0 

c 

rxmin 

= 

rkl*dexp (-thl/tkO) 

c 

qdetO 

= 

4 . 530d5 

c 

RDX  Test 

pcj 

= 

2  6 . 5d9 

rkl 

= 

110d6 

rk2 

= 

OdO 

pexp 

= 

3 . 5d0 

zexp 

= 

0 . 93d0 

thl 

= 

OdO 

th2 

= 

OdO 

rxmin 

= 

rkl* ( (pamb/pcj ) **pexp) 

qdetO 

— 

5 . 375d6 

c 

Particle 

data 

xpl 

= 

1.0d-2 

xp2 

= 

5 . 9d-2 

pmass 

= 

4 . 3d0 

rop 

= 

7  8  60d0 

rdp 

= 

280d-6 

pep 

= 

446d0 

mu 

= 

1 . 7d-5 

c 

mu 

= 

1.0d-3 

tcon 

= 

2 . 57d-2 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Code  control  data  and  flags 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Data 


off 

=  ld-6 

cfl 

=  0 . 5d0 

n 

=  0 

nfil 

=  0 

nstart 

=  0 

nstp 

=  10 

ndmp 

=  1 

dtmx 

=  ld-2 

time 

=  OdO 

tend 

=  50d0 

Flags 

irst 

=  1 

iav 

=  1 

iext 

=  1 

ilim 

=  1 

ieos 

=  3 

igeo 

=  1 

irxn 

=  1 

iefx 

=  2 

ipar 

=  0 

idrg 

=  1 

imach 

=  1 
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c  Debug  control 
idbgl  =  0 
idbgf  =  0 
idbgs  =  0 
idbgp  =  0 


write  (  * , * 

i 

Code 

Control  Data 

write  (  * , * 

i 

nstp 

= 

' , nstp 

write  (  * , * 

i 

ndmp 

= 

' , ndmp 

write  (  * , * 

i 

tend 

= 

' , tend 

if  (ipar 

eq 

1 )  write ( * , * )  ' 

write  (  * , * 

1 

1 

write  (  * , * 

1 

Flags 

.  I 

write  (  * , * 

1 

irst 

= 

' , irst 

write  (  * , * 

1 

iav 

= 

'  ,  iav 

write  (  * , * 

1 

iext 

= 

' , iext 

write  (  * , * 

1 

ilim 

= 

' , ilim 

write  (  * , * 

1 

1 

write  (  * , * 

1 

ieos 

= 

' , ieos 

write  (  * , * 

I 

igeo 

= 

' , igeo 

write  (  * , * 

1 

irxn 

= 

' , irxn 

write  (  * , * 

1 

iefx 

= 

' , iefx 

write  (  * , * 

1 

i 

write  (  * , * 

1 

ipar 

= 

' , ipar 

write  (  * , * 

1 

idrg 

= 

' , idrg 

write  (  * , * 

1 

imach 

= 

' , imach 

write  (  * , * 

1 

1 

pause 

c  Derived  data 

c  Thermal  data 

cpg  =  rgas 

+  cvg 

ppr  =  cpg*tcon/mu 
crppr  =  ppr**cl3 


c  EOS  Parameters 


rhl 

=  rl*r0 

rh2 

=  r2*r0 

wrl 

=  wj/rhl 

wrlr 

=  wrl/rO 

wr2 

=  w j / rh2 

wr2r 

=  wr2/r0 

cjh 

=  cj  * (rO** 

nhpl 

=  nh  +  ldO 

alf  a 

=  nh  -  ldO 

nhml 

=  alf  a 

nhm2 

=  nh  -  2d0 

eO 

=  cvg*tk0 

'  ,  npar 


c  Hayes-I  EOS 

t3  =  cvs*tkO*gh/rO 

t4  =  hi /rO /nh/alf a 

t5  =  pamb/gh  +  t4 

t7  =  pamb/gh  +  beta  +  t4 
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thta  =  t3  -  pamb/rO 
beta  =  thta  +  alfa*t4 


c  Compute  coefficients  for  Hayes  pressure  derivatives 


C  (1) 

=  ldO 

c  (2 ) 

=  beta 

c  (3) 

=  -t4 

c  (4) 

=  t5 

c  (5) 

=  aj 

c  (6) 

=  bj 

c  (7) 

=  qdetO  +  eO 

c  (8) 

=  hl/gh/nh 

c 


c 


Particle  phase  parameters 
dip  =  2d0*rdp 

pOmas  =  c43*pi*rop*rdp*rdp*rdp 
pvol  =  pmass/rop 
if  (chx  .le.  x2)  then 
write ( * , * )  '  ' 

write (*,*)  '  chx  <  x2 . ' 

write ( *  ,  *  )  '  ' 

stop 
else 

dx  =  chx  -  x2 
endif 

cvol  =  c43*pi*x2*x2*x2 
cvol  =  hevol  +  pvol 
alf2  =  pvol/cvol 
alf 1  =  ldO  -  alf 2 


if  (ipar  .eq.  1  .and.  alfl  .eq.  OdO)  then 
write ( * , * )  '  ' 

write (*,*)  '  alfl  =0!' 

write ( * , * )  '  ' 

stop 
endif 

alf 21  =  alf 2/alf 1 

c  Other  constants 

epsm  =  cl4*(ld0  -  kap) 
epsp  =  cl4* (ldO  +  kap) 
garni  =  gamm  -  ldO 

c  Set  up  the  solver  report  file 

open ( 90 , f ile= ' rpt . txt ' , f orm= ' formatted ' ) 

write  (90,*)  '  **********  Detonation  Solver  Report  File 

•k'k'k'k'k'k'k'k'k'k  * 


write ( 90 , * ) 

1  1 

write  (  90 , * ) 

'  Reaction 

Data :  ' 

write ( 90 , * ) 

'  qdet  = 

' , qdetO 

write  (  90 , * ) 

'  eact  = 

'  ,  eact 

write ( 90 , * ) 

'  rk  = 

'  ,  rk 

write  (  90 , * ) 

'  rkl  = 

'  ,  rkl 

write ( 90 , * ) 

'  rk2  = 

'  ,  rk2 

write  (  90 , * ) 

'  pexp  = 

'  ,  pexp 

write ( 90 , * ) 

'  zexp  = 

' , zexp 
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write ( 90 , * ) 

'  Pc  j  = 

'  ,pcj 

write  (  90 , * ) 

'  thl 

'  ,  thl 

write ( 90 , * ) 

'  th2 

'  ,  th2 

write  (  90 , * ) 

I  1 

write ( 90 , * ) 

'  rxmin  = 

' , rxmin 

write  (  90 , * ) 

i  i 

write ( 90 , * ) 

'  EOS  Control  Data 

write  (  90 , * ) 

'  ztoll  = 

' , ztoll 

write ( 90 , * ) 

'  ztol2  = 

' , ztol2 

write  (  90 , * ) 

1  I 

write ( 90 , * ) 

'  CPG  EOS 

Data :  ' 

write  (  90 , * ) 

'  gamm  =  ' 

,  gamm 

write ( 90 , * ) 

'  gaml  =  ' 

,  gaml 

write  (  90 , * ) 

1  1 

write ( 90 , * ) 

'  Hayes-I 

EOS  Data 

write  (  90 , * ) 

'  HI 

'  r  hi 

write ( 90 , * ) 

'  Cvs  = 

'  ,  CVS 

write  (  90 , * ) 

'  g 

-,gh 

write ( 90 , * ) 

'  N 

'  ,nh 

write  (  90 , * ) 

'  TO 

'  ,  tkO 

write ( 90 , * ) 
do  nn  =  1,8 

1  1 

write  (  90 , * 
enddo 

)  '  c ( ' , nn 

'  )  =  ' 
r  )  r 

write ( 90 , * ) 

1  1 

write  (  90 , * ) 

'  alf a  =  ' 

,  alf  a 

write ( 90 , * ) 

'  beta  =  ' 

,  beta 

write  (  90 , * ) 

'  thta  =  ' 

,  thta 

write ( 90 , * ) 

'  t3 

,  t3 

write  (  90 , * ) 

'  t4 

,  t4 

write ( 90 , * ) 

'  t5 

,  t5 

write  (  90 , * ) 

'  t7 

,t7 

write ( 90 , * ) 

1  I 

write  (  90 , * ) 

'  JWL  EOS 

Data :  ' 

write ( 90 , * ) 

'  rO 

,  rO 

write  (  90 , * ) 

'A  =  ' 

,  aj 

write ( 90 , * ) 

'  B  =  ' 

,bj 

write  (  90 , * ) 

'  C  =  ' 

,  cj 

write ( 90 , * ) 

'  R1 

,  rl 

write  (  90 , * ) 

'  R2 

,  r2 

write ( 90 , * ) 

'  W  =  ' 

,  wj 

write  (  90 , * ) 

'  Cvg  =  ' 

,  cvg 

write ( 90 , * ) 

'  Cpg  =  1 

,  cpg 

write  (  90 , * ) 

'  eO 

,  eO 

write ( 90 , * ) 

1  1 

write  (  90 , * ) 

'  Particle 

Data :  ' 

write ( 90 , * ) 

'  pmass  = 

'  ,  pmass 

write  (  90 , * ) 

'  rop  = 

'  ,  rop 

write ( 90 , * ) 

'  rdp  = 

'  ,  rdp 

write  (  90 , * ) 

'  dip 

'  ,  dip 

write ( 90 , * ) 

'  mu  = 

'  ,  mu 

write  (  90 , * ) 

'  tcon  = 

' , tcon 

write ( 90 , * ) 

'  ppr  = 

'  i  ppr 

write  (  90 , * ) 

'  pOmas  = 

'  ,  pOmas 

write ( 90 , * ) 

'  hevol  = 

'  ,  hevol 

write  (  90 , * ) 

'  pvol  = 

'  ,  pvol 

write ( 90 , * ) 

'  cvol  = 

'  ,  cvol 
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alf  1 
alf  2 


alf  1 
alf  2 


write  (  90 , * ) 
write ( 90 , * ) 
write  (  90 , * ) 
write  (90,*)  '  Other  Data:' 

write (90,*)  '  kap  =  ' , kap 

write (90,*)  '  epsm  =  ' , epsm 

write  (90,*)  '  epsp  =  ' , epsp 

write  (  90 , * )  '  ' 

close  ( 90 ) 

write  ( * , * )  '  ' 

write  (*,*)  '  Report  file  ready.' 

write  (  * , * )  '  ' 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Grid  Generation  Section 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
dx  =  (x2  -  xl) / (imax-1) 
offs  =  O.ldO 
do  i  =  1 , imax 

x(i)  =  xl  +  (i-l)*dx 

c  write(*,*)  '  i  =  ' , i , '  x  =  ',x(i) 

enddo 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Initial  Conditions  and  Restart  File  Section 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Set  initial  conditions  (no  restart) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
if  (irst  .eq.  0)  then 

c  Set  time  zero  primitive  variables 
do  i  =  1 , imax-1 

xc  =  cl2*(x(i)  +  x(i+l)) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  CPG  EOS  ICs 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
if  (ieos  .eq.  0)  then 

r(i)  =  Id0/(ld0  +  3d0*dexp (-xc*xc) ) 

p(i)  =  ldO 

u ( i )  =  OdO 

z ( i )  =  OdO 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  JWL  EOS  ICs 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
else  if  (ieos  .eq.  1)  then 
r ( i )  =  1 . 2d0 

c  p(i)  =  25d0*pamb/ (1 . OOOOldO  -  dexp (-xc*xc) ) 
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P  ( i ) 


=  ( (x2-xc) * (40d0*pamb/ (1 . OOOOldO  -  dexp (-xc*xc) ) ) 
&  +  x2*pamb)/x2 

c  write(70,*)  xc, '  ',p(i) 

c  if  (xc  .It.  offs)  then 

c  p(i)  =  fct* (xc-offs) * (xc-offs)  +  pamb 

c  else 

c  p ( i )  =  pamb 

c  endif 

u ( i )  =  OdO 

z ( i )  =  OdO 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Hayes-I/JWL  EOS  ICs 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
else  if  (ieos  .eq.  2)  then 
r  ( i )  =  rO 

c  p(i)  =  25d0*pamb/ (1 . OOOOldO  -  dexp (-xc*xc) ) 

c  p(i)  =  ( (x2-xc) * (25d0*pamb/ (1 . OOOOldO  -  dexp (-xc*xc) ) ) 

c  &  +  x2*pamb)/x2 

c  p ( i )  =  pamb 

c  HMX  or  NM 

p(i)  =  2d0*pc j *dexp ( -xc*xc/0 . 00 ld0/0 . 00 ldO )  +  pamb 

u ( i )  =  OdO 

z ( i )  =  OdO 

tk(i)  =  tkO 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Hayes-II/JWL  EOS  ICs 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
else  if  (ieos  .eq.  3)  then 
r  ( i )  =  rO 

c  HMX/RDX/NM 

c  p(i)  =  2d0*pcj *dexp (-xc*xc/0 . 004d0/0 . 004d0)  +  pamb 

if  (i  .le.  100)  then 

p(i)  =  5d0*pcj  +  pamb 

else 

p ( i )  =  pamb 

endif 

c  NM 

c  p(i)  =  2d0*pcj *dexp (-xc*xc/0 . 0005d0/0 . 0005d0)  +  pamb 

u ( i )  =  OdO 

z ( i )  =  OdO 

tk(i)  =  (p(i)  -  pamb) /cvs/gh  +  tkO 
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else 

write  ( * , * )  '  ' 

write  (*,*)  '  Unknown  EOS' 

write  ( * , * )  '  ' 

stop 
endif 

enddo 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Particle  ICs 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
if  (ipar  .eq.  1)  then 

c  Check  particle  bounds 

if  (xpl  .It.  xl  .or.  xp2  . gt.  x2)  then 
write  ( * , * )  '  ' 

write  (*,*)  '  Particle  X  limits  are  wrong.' 

write  ( * , * )  '  ' 

stop 
endif 

dxp  =  (xp2  -  xpl)/ (npar  -  1) 
do  np  =  l,npar 

px(np)  =  xpl  +  (np-l)*dxp 
pu (np)  =  OdO 
ptk(np)  =  tkO 
pa (np)  =  OdO 
pq(np)  =  OdO 

write ( * , * )  px(np),'  ’,pu(np),'  ’,pa(np) 
enddo 
pause 

write  ( * , * )  '  ' 

write (*,*)  '  Particles  ready.' 

write  (  * , * )  '  ' 

endif 

else  if  (irst  .eq.  1)  then 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Read  the  restart  file 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
write (*,*)  '  Reading  restart  file.' 

open ( 40 , f ile= ' restart . data ' , f orm= ' unformatted ' ) 
read(40)  nstart 
read(40)  nfil 
read(40)  time 
do  i  =  1 , imax-1 

read (40)  r ( i ) , p ( i ) , u ( i ) , z  ( i ) 
enddo 
close  (40) 

else 

write  ( * , * )  '  ' 

write  (*,*)  '  Unknown  restart  option.' 
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write  (  * , * )  '  ' 

endif 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Compute  initial  derived  flow  variables  for  the  cells 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
do  i  =  1 , imax-1 

if  (ieos  .eq.  0)  then 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  CPG  EOS  internal  energy  and  pressure  derivatives 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ei(i)  =  p (i) /r (i) /gaml  -  z(i)*qdet0 

dpdr  =  gaml*ei(i)  +  gaml*z (i) *qdet0 

dpde  =  gaml*r(i) 

dpdz  =  gaml*r (i) *qdet0 

else  if  (ieos  .eq.  1)  then 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  JWL  EOS  internal  energy  and  pressure  derivatives 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
rht  =  r ( i ) /r0 
rhti  =  ldO/rht 
ri  =  ldO /r ( i ) 

tmp  =  p(i)  -  aj*(ld0  -  wrl*r (i) ) *dexp (-rhl*ri) 

&  -  bj*(ld0  -  wr2*r (i) ) *dexp (-rh2*ri) 

ei(i)  =  tmp/wj*ri  -  z(i)*qdet0 

tmp  =  a j * (rhl*ri*ri  -  wj*ri  -  wj /rhl ) *dexp ( -rhl*ri ) 

tmp  =  tmp  +  bj * (rh2*ri*ri  -  wj*ri  -  wj /rh2 ) *dexp ( -rh2*ri ) 

dpdr  =  tmp  +  wj*ei(i)  +  wj * z ( i ) *qdet0 

dpde  =  wj*r(i) 

dpdz  =  wj *r ( i ) *qdet0 

else  if  (ieos  .eq.  2)  then 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Hayes-I/JWL  EOS  internal  energy  and  pressure  derivatives 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


ra 

=  r  (i) 

ra2 

=  ra*ra 

za 

=  z  (i) 

rz 

=  ra*za 

omz 

=  ldO  - 

c  Solid  phase  limit 

if  (za  .le.  ztoll)  then 

ei(i)  =  p(i)/gh  +  beta*r0/ra  +  t4* ( (ra/rO) **alfa)  -  t7 
dpdr  =  beta*r0*gh/ra2 
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& 


-  alfa*gh*t4* (ra** (alfa-ldO) ) / (rO**alfa) 


dpde  =  gh 
c  Mixed  phases 

else  if  (ztoll  .It.  za  .and.  za  .It.  ztol2)  then 


c  Evaluate  denominator  functions 

bot  =  omz/gh  +  ldO/wj/ra 
if  (bot  .It.  Id-10)  then 
write  (  * , * )  '  ' 

write (*,*)  '  Zero  denonimator  term.' 

write  (  * , * )  '  ' 

stop 
endif 

bot2  =  bot*bot 
botr  =  -Id0/wj/ra2 


c  Evaluate 


numerator  functions 


top 

(2) 

=  omz  -  rO/ra 

top 

(3) 

=  (omz* *nh) *( (ra/rO ) ** 

alfa) 

top 

(4) 

=  omz 

top 

(5) 

=  (ldO/wj/ra  -  za/rhl) 

*dexp ( -rhl / rz 

top 

(6) 

=  (ldO/wj/ra  -  za/rh2) 

*dexp ( -rh2 / rz 

top 

(7) 

=  za 

c  Compute  internal  energy 

ei ( i )  =  bot*p ( i ) 
do  nn  =  2,7 

ei(i)  =  ei(i)  -  c (nn) *top (nn) 
enddo 

top ( 1 )  =  ei  ( i ) 


c  Compute  derivatives  for  numerator  functions 


topr ( 1 ) 

=  OdO 

topr (2 ) 

=  r0/ra2 

topr  ( 3 ) 

=  alfa/rO* (omz* *nh) *( (ra/rO ) 

** (alfa-ldO ) ) 

topr ( 4 ) 

=  OdO 

topr  ( 5 ) 

=  (rhl/wj/rz  -  ldO/wj 

-  ldO ) 

*dexp ( -rhl/rz ) 

/ra2 

topr  ( 6 ) 

=  (rh2/wj/rz  -  ldO/wj 

-  ldO ) 

*dexp ( -rh2 /rz ) 

/ra2 

topr ( 7 ) 

=  OdO 

c  Compute  density  and  internal  energy  derivatives  of  pressure 
dpdr  =  OdO 
do  nn  =  1,7 

dpdr  =  dpdr  +  c (nn) * (bot*topr (nn)  -  botr*top (nn) ) 
enddo 

dpdr  =  dpdr/bot2 
dpde  =  ldO/bot 

c  Gas  phase  limit 
else 

ei  (i)  =  p  (i)  /wj  /ra 

&  -  a j * ( ldO /wj /ra  -  ldO /rhl ) *dexp ( -rhl /ra) 

&  -  b j * ( ldO /wj /ra  -  Id0/rh2 ) *dexp ( -rh2/ra) 

&  -  qdetO  -  eO 
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& 

& 

& 


dpdr 


=  wj  *ei ( i ) 

+  aj*(rhl/ra2  -  w j /ra  -  wj /rhl ) *dexp ( -rhl /ra) 
+  bj*(rh2/ra2  -  w j /ra  -  wj /rh2 ) *dexp ( -rh2/ra) 
+  w j  * (qdetO  +  eO ) 

dpdz  =  aj*(rhl/ra  -  wj  -  wj *ra/rhl ) *dexp ( -rhl /ra) 

&  +  bj*(rh2/ra  -  wj  -  wj *ra/rh2 ) *dexp ( -rh2/ra) 

&  +  ra*wj* (qdetO  +  eO) 

dpde  =  wj*ra 

endif 

else  if  (ieos  .eq.  3)  then 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Hayes-II/JWL  EOS  internal  energy  and  pressure  derivatives 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


ra 

=  r  (i) 

ra2 

=  ra*ra 

za 

=  z  (i) 

rz 

=  ra*za 

omz 

=  ldO  - 

c  Solid  phase  limit 

if  (za  .le.  ztoll)  then 


& 

ei  (i) 

=  p(i)/gh  +  beta*r0/ra  +  t4* 
-  hi /gh/nh* (( (ra/rO ) **nh)  - 

( (ra/rO ) **alfa)  -  t7 
ldO ) 

& 

& 

dpdr 

=  beta*r0*gh/ra2 
-  alfa*gh*t4* (ra** (alfa-ldO) 
+  hi /r0* ( (ra/rO ) **nhml ) 

) / (r0**alfa) 

dpde 

=  gh 

Mixed 

phases 
else  if 

(ztoll  .It.  za  .and.  za  .It. 

ztol2)  then 

c  Evaluate  denominator  functions 

bot  =  omz/gh  +  IdO/wj/ra 
if  (bot  .it.  Id-10)  then 
write  (  * , * )  '  ' 

write (*,*)  '  Zero  denonimator  term.' 

write  (  * , * )  '  ' 

stop 
endif 

bot2  =  bot*bot 
botr  =  -Id0/wj/ra2 

c  Evaluate  numerator  functions 

top (2)  =  omz  -  rO/ra 

top(3)  =  (omz**nh) * ( (ra/rO ) **alfa) 

top (4)  =  omz 

top(5)  =  (IdO/wj/ra  -  za/rhl ) *dexp ( -rhl /rz ) 

top(6)  =  (IdO/wj/ra  -  za/rh2 ) *dexp ( -rh2 /rz ) 
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top (7) 
top ( 8 ) 


=  za 

=  (omz**nhpl )*( (ra/rO ) **nh)  +  za  -  ldO 

c  Compute  internal  energy 

ei ( i )  =  bot*p ( i ) 
do  nn  =  2,8 

ei(i)  =  ei(i)  -  c (nn) *top (nn) 
enddo 

top ( 1 )  =  ei  ( i ) 

c  Compute  derivatives  for  numerator  functions 


topr 

(1) 

=  OdO 

topr 

(2) 

=  r0/ra2 

topr 

(3) 

=  alf a/rO* (omz* *nh) *(( ra/rO  ) 

** (alfa-ldO ) ) 

topr 

(4) 

=  OdO 

topr 

(5) 

=  (rhl/wj/rz  -  ldO/wj  - 

-  ldO ) 

*dexp ( -rhl /rz ) 

/ra2 

topr 

(6) 

=  (rh2/wj/rz  -  ldO/wj  - 

-  ldO ) 

*dexp (-rh2/ rz) 

/ra2 

topr 

(7) 

=  OdO 

topr 

(8) 

=  nh/rO* (omz**nhpl )*( (ra/rO ) 

**nhml ) 

c  Compute  density  and  internal  energy  derivatives  of  pressure 
dpdr  =  OdO 
do  nn  =  1,8 

dpdr  =  dpdr  +  c (nn) * (bot*topr (nn)  -  botr*top (nn) ) 
enddo 

dpdr  =  dpdr/bot2 
dpde  =  IdO/bot 

c  Gas  phase  limit 
else 

ei  (i)  =  p  (i)  /wj  /ra 

&  -  a j * ( ldO /wj /ra  -  ldO/rhl ) *dexp ( -rhl /ra) 

&  -  b j * ( ldO /wj /ra  -  Id0/rh2 ) *dexp ( -rh2/ra) 

&  -  qdetO  -  eO 

dpdr  =  wj  *ei ( i ) 

&  +  aj*(rhl/ra2  -  w j /ra  -  wj /rhl ) *dexp ( -rhl /ra) 

&  +  bj*(rh2/ra2  -  w j /ra  -  wj /rh2 ) *dexp ( -rh2/ra) 

&  +  wj* (qdetO  +  eO) 

dpdz  =  aj*(rhl/ra  -  wj  -  wj *ra/rhl ) *dexp ( -rhl /ra) 

&  +  bj*(rh2/ra  -  wj  -  wj *ra/rh2 ) *dexp ( -rh2/ra) 

&  +  ra*wj* (qdetO  +  eO) 

dpde  =  wj*ra 

endif 

else 

write  ( * , * )  '  ' 

write  (*,*)  '  Unknown  EOS' 

write  (  * , * )  '  ' 

stop 
endif 

c  Compute  the  speed  of  sound 
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if  (dpdr  .It.  OdO)  dpdr  =  dabs (dpdr) 
a2  =  dpdr  +  p (i) *dpde/r (i) /r (i) 

if  (a2  .It.  OdO)  then 
write  (  * , * )  '  ' 

write  (*,*)  '  Negative  initial  squared  sound  speed!' 

write  (  * , * )  '  i  =  ' , i 

write  (  * , * )  '  ' 

stop 
endif 

a (i)  =  dsqrt (a2) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Initial  reaction  rate 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Floor  on  1  -  z  near  0 

if  (z(i)  .gt.  ztol2)  then 
omz  =  OdO 
else 

omz  =  IdO  -  z ( i ) 
endif 

c  Test  Rate  1 

c  rxr(i)  =  rkl*dsqrt (omz) 

c  if  (p(i,j)  -  ld9  .It.  OdO)  rxr(i)  =  OdO 

c  if  (p(i,j)  -  ld9  .eq.  OdO)  rxr(i)  =  0.5d0*rxr(i) 

c  CPG  Test  Rate 

c  rxr(i)  =  rk*omz*dexp (-eact*r (i) /p (i) )  -  rxmin 

c  HMX  Test  Rate 

c  rxr(i)  =  rkl* (omz**zexp) * ( (p (i) /pcj ) **pexp)  -  rxmin 

c  if  (rxr(i)  .It.  OdO)  rxr(i)  =  OdO 

c  NM  Test  Rate 

c  rxr(i)  =  (rkl*dexp (-thl/tk (i) ) *omz 

c  &  +  rk2*dexp (-th2/tk (i) ) *z (i) ) * (omz**zexp)  -  rxmin 

c  if  (rxr(i)  .It.  OdO)  rxr(i)  =  OdO 

c  RDX  Test  Rate 

rxr(i)  =  rkl* (omz**zexp) * ( (p (i) /pcj ) **pexp)  -  rxmin 
if  (rxr(i)  .It.  OdO)  rxr(i)  =  OdO 

enddo 

c  Write  the  initial  conditions  files 
if  (irst  .eq.  0)  then 

open (21 , f ile= ' heic . dat ' , form= ' formatted ' ) 

7  0  format ( lx,  dl2 . 6 ,  lx,dl2.6, lx,dl2.6, lx,dl2.6, lx,dl2.6, lx,dl2.6, 

&  lx,  dl2 . 6,  lx,  dl2 . 6) 

72  format (Ix,dl2.6,lx,dl2.6,lx,dl2.6,lx,dl2.6,lx,dl2.6,lx,dl2.6, 

&  lx,  dl2 . 6,  lx,  dl2 . 6,  lx,  dl2 . 6) 

do  i  =  1 , imax-1 

xc  =  cl2*(x(i)  +  x(i+l)) 

write (21,72)  xc,  r (i)  ,  u (i)  ,  p (i) , z (i) , ei (i) , a (i) , rxr (i) , tk (i) 
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enddo 
close (21) 

write  (*,*)  '  ICs  ready.' 

write  ( * , * )  '  ' 

if  (ipar  .eq.  1)  then 

open (21 , f ile= ' paic . dat ' , form= ' formatted ' ) 
do  np  =  l,npar 

write(21,*)  px(np),'  ',0d0,'  ' ,pu(np) 

enddo 
close (21) 
endif 
endif 

pause 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Set  the  internal  energy  correction  and  scale  variables 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
eOcr  =  OdO 
eta  =  0.999d0 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Main  Solver  Loop 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
do  while  (n  .It.  nstp  .and.  time  .It.  tend) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Allocate  particles  to  cells 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
if  (ipar  .eq.  1)  then 
do  np  =  l,npar 


peel (np)  = 

int((px(np)  -  xl)/dx)  +  1 

c 

write ( * , * ) 

'  px ( ' , np, ' )  =  ' , px (np) 

c 

write ( * , * ) 

'  pcel(',np, ')  =  ' ,pcel(np) 

c 

write ( * , * ) 

1  1 

enddo 


pum  =  OdO 
pam  =  OdO 
do  np  =  l,npar 

if  (dabs (pu (np) )  . gt.  pum)  pum  =  dabs (pu (np) ) 
if  (dabs (pa (np) )  .gt.  pam)  pam  =  dabs (pa (np) ) 
enddo 
endif 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Compute  time  step 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
dt  =  ld2 
do  i  =  1 , imax-1 

dx  =  x(i+l)  -  x ( i ) 

dtO  =  dx/ (dabs (u  ( i ) )  +  a(i)) 
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if  (ipar  .eq.  1)  then 

dtl  =  dx/ (dabs (u ( i ) )  +  pum) 

dtO  =  min (dtO , dtl ) 
c  dtl  =  2dl*dx/pam 

c  dtO  =  min (dtO , dtl ) 

endif 

if  (dtO  .It.  dt)  dt  =  dtO 
enddo 

dt  =  cfl*dt 

dt  =  min(dt,dtmx) 

if  (idbgl  .eq.  1)  then 
write (*,*)  '  dt  =  ' , dt 

write ( * , * )  '  ' 

endif 

c  Set  boundary  conditions 


c 

Symmetric  at 

x  =  0 

r  (0) 

r  (1) 

u  ( 0  )  = 

— u  ( 1 ) 

p  (0)  = 

P  (1) 

z  (0) 

Z  (1) 

ei  ( 0 )  = 

ei  (1) 

c 

Fixed  at  x  = 

xmax 

c 

r ( imax) 

=  ldO 

c 

u ( imax) 

=  OdO 

c 

p ( imax) 

=  ldO 

c 

z ( imax) 

=  OdO 

c 

ei ( imax 

)  =  p ( imax) /r ( imax) /garni 

c 

Extrapolated 

at  x  =  xmax 

r ( imax) 

=  r (imax-1) 

u ( imax) 

=  u (imax-1) 

p ( imax) 

=  p (imax-1) 

z ( imax) 

=  z (imax-1) 

ei(imax)  =  ei(imax-l) 

if  (idbgl  .eq.  1)  then 
write  (  * , * )  '  BCs :  ' 

write(*,*)  '  r(0)  =  ',r(0) 

write(*,*)  '  u(0)  =  ',u(0) 

write(*,*)  '  p(0)  =  ',p(0) 

write(*,*)  '  z(0)  =  ',z(0) 

write(*,*)  '  ei(0)  =  ',ei(0) 

write ( * , * )  '  ' 

write (*,*)  '  r(imax)  =  ' ,r(imax) 

write (*,*)  '  u(imax)  =  ' ,u (imax) 

write(*,*)  '  p(imax)  =  ' ,p(imax) 
write  (*,*)  '  z(imax)  =  1 , z (imax) 

write (*,*)  '  ei(imax)  =  ',ei(imax) 
write ( * , * )  '  ' 

endif 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Minimum  reaction  rate  taken  at  cell  imax-1 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
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c  Floor  on  1  -  z  near  0 

if  (z(imax-l)  .gt.  ztol2)  then 
omz  =  OdO 
else 

omz  =  IdO  -  z(imax-l) 
endif 

c  HMX  or  RDX  Test 

rxmin  =  rkl* (omz**zexp) *( (p (imax-1) /pcj ) **pexp) 
c  NM  Test 

c  rxmin  =  (rkl*dexp (-thl/tk (imax-1) ) *omz 

c  &  +  rk2*dexp (-th2/tk (imax-1) ) *z (imax-1) )* (omz**zexp) 

c  write (*,*)  '  rxmin  =  rxmin 

c  write (*,*)  '  rxr  =  rxr ( imax-1 ) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Compute  conservative  variables;  assemble  source  terms 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
do  i  =  1 , imax-1 


qv ( i , 1 ) 

=  r  (i) 

qv ( i , 2 ) 

=  r  (i)  *u  (i) 

et 

=  ei(i)  +  0.5d0*u(i 

qv ( i , 3 ) 

=  r ( i ) *et 

qv ( i , 4 ) 

=  r  (i)  *z  (i) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Compute  the  source  vectors 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Geometric 

xc  =  cl2*(x(i)  +  x(i+l)) 

sg(i,l)  =  -r ( i ) *u ( i ) /xc 
sg(i,2)  =  -r ( i ) *u ( i ) *u ( i ) /xc 

sg(i,3)  =  -u (i) * (r  (i) *et  +  p(i))/xc 

sg(i,4)  =  -r ( i ) *u ( i ) * z ( i ) /xc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Reaction  rate 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Floor  on  1  -  z  near  0 

if  (z(i)  .gt.  ztol2)  then 
omz  =  OdO 
else 

omz  =  IdO  -  z ( i ) 
endif 

c  CPG  Test  Rate 

c  rxr(i)  =  rk*omz*dexp (-eact*r (i) /p (i) )  -  rxmin 

c  HMX  Test  Rate 

c  rxr(i)  =  rkl* (omz**zexp) * ( (p (i) /pcj ) **pexp)  -  rxmin 

c  if  (rxr(i)  .It.  OdO)  rxr(i)  =  OdO 

c  NM  Test  Rate 
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c 

rxr ( i )  = 

( rkl *dexp ( -thl /tk ( i )  ) 

)  *omz 

c 

& 

+ 

rk2 *dexp ( -th2 /tk ( i )  ] 

1  *z (i) ) * (omz**zexp)  -  rxmin 

c 

if  (rxr  (i 

)  .It.  OdO)  rxr(i)  = 

OdO 

c  RDX  Test  Rate 

rxr(i)  =  rkl* (omz**zexp) * ( (p (i) /pcj ) **pexp)  -  rxmin 
if  (rxr(i)  .It.  OdO)  rxr(i)  =  OdO 


c  Reaction  rate  terms 


srx ( i , 1 ) 

=  OdO 

srx ( i , 2 ) 

=  OdO 

srx ( i , 3 ) 

=  OdO 

srx ( i , 4 ) 

=  r  ( i 

phase 

sp (i, 1) 

=  OdO 

sp (i, 2) 

=  OdO 

sp  (i,  3) 

=  OdO 

sp (i, 4) 

=  OdO 

enddo 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Compute  particle  phase  coupling  terms 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
if  (ipar  .eq.  1)  then 
do  np  =  l,npar 


c  Momentum 

sp (peel (np) , 2) 

c  Energy 

sp (peel (np) , 3) 


sp (peel (np) , 2 )  -  pOmas*pa (np) 

sp (peel (np) , 3)  -  pq(np) 


enddo 

endif 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Compute  the  total  source  vector 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
do  i  =  1 , imax-1 

c  if  (sp(i,2)  .ne.  OdO)  then 

c  write(*,*)  '  i  =  ' , i , '  sp  =  ',sp(i,2) 

c  endif 

do  m  =  1,4 

s(i,m)  =  igeo*sg(i,m)  +  irxn*srx (i,m)  +  ipar*sp(i,m) 


enddo 

if  (idbgs  .eq. 

1)  then 

write ( * , * )  ' 

i  =  '  ,  i 

write ( * , * )  ' 

ql  =  '  ,  qv  (i,  1 

write  (  * , * )  ' 

q2  =  ' , qv ( i , 2 

write  (  * , * )  ' 

q3  =  ' , qv (i, 3 

write  (  * , * )  ' 

q4  =  ' , qv ( i , 4 

write ( * , * )  ' 

1 
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write ( * , * ) 

'  si  = 

'  ,  s (i,  1 

write ( * , * ) 

'  s2  = 

'  ,  s (i, 2 

write  (  * , * ) 

'  s3  = 

'  ,  s  ( i ,  3 

write  (  * , * ) 

II 

U) 

'  ,  s (i,  4 

write ( * , * ) 

i  i 

pause 

endif 

enddo 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Compute  the  numerical  flux  at  each  grid  point 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
71  format (2x,dl2.6,2x,dl2.6,2x,dl2.6,2x,dl2.6) 

do  i  =  1 , imax 

c  Left  interface  variables 

if  (i  .eq.  1)  then 
c  First  order  at  the  boundary 
rl  =  r(i-l) 
ul  =  u ( i — 1 ) 
pi  =  p ( i — 1 ) 
zl  =  z  (i-1) 

rr  =  r ( i ) 
ur  =  u  ( i ) 
pr  =  p (i) 
zr  =  z (i) 

else  if  (2  .le.  i  .and.  i  .le.  imax-1)  then 
c  Higher-order 

if  (ilim  .eq.  0)  then 

c  Hossaini  limiting  strategy 

dqwr  =  r(i-l)  -  r(i-2) 

dqer  =  r(i)  -  r(i-l) 

dqir  =  r(i+l)  -  r(i) 

phir  =  cl4* (2d0*dqwr*dqer  +  eps) 

&  / (dqwr*dqwr  +  dqer*dqer  +  eps) 

dqwu  =  u(i-l)  -  u(i-2) 

dqeu  =  u(i)  -  u(i-l) 

dqiu  =  u(i+l)  -  u(i) 

phiu  =  cl4* (2d0*dqwu*dqeu  +  eps) 

&  / (dqwu*dqwu  +  dqeu*dqeu  +  eps) 

dqwp  =  p (i-1)  -  p  ( i-2 ) 

dqep  =  p ( i )  -  p(i-l) 

dqip  =  p  ( i  +  1 )  -  p ( i ) 

phip  =  cl4* (2d0*dqwp*dqep  +  eps) 

&  / (dqwp*dqwp  +  dqep*dqep  +  eps) 

dqwz  =  z(i-l)  -  z(i-2) 

dqez  =  z(i)  -  z(i-l) 

dqiz  =  z(i+l)  -  z(i) 

phiz  =  cl4* (2d0*dqwz*dqez  +  eps) 
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& 


/ (dqwz*dqwz  +  dqez*dqez  +  eps) 


c  Density 


c  Velocity 


c  Pressure 


c  Rx  Progress 


rl 

=  r  ( i  —  1 ) 

+ 

iext*phir* 

(epsm*dqwr 

+ 

epsp*dqer ) 

rr 

=  r  (i) 

iext*phir* 

(epsm*dqir 

+ 

epsp*dqer ) 

ul 

=  u(i-l) 

+ 

iext*phiu* 

(epsm*dqwu 

+ 

epsp*dqeu) 

ur 

=  u(i) 

iext*phiu* 

(epsm*dqiu 

+ 

epsp*dqeu) 

Pi 

=  P(i-l) 

+ 

iext*phip* 

(epsm*dqwp 

+ 

epsp*dqep) 

pr 

=  P  ( i ) 

iext*phip* 

(epsm*dqip 

+ 

epsp*dqep) 

zl 

=  z  ( i  —  1 ) 

+ 

iext*phiz* 

(epsm*dqwz 

+ 

epsp*dqez ) 

zr 

=  z  (i) 

- 

iext*phiz* 

(epsm*dqiz 

+ 

epsp*dqez ) 

else  if  (ilim  .eq.  1)  then 
c  Hirsch  limiting  strategy 


dra 

= 

r  ( i  +  1 ) 

- 

r  (i) 

drb 

= 

r  (i) 

- 

r  ( i-1 

drc 

= 

r  ( i  —  1 ) 

- 

r  ( i-2 

drd 

= 

drb 

- 

drc 

dre 

= 

dra 

- 

drb 

dua 

= 

u  (i+1) 

_ 

u  ( i ) 

dub 

= 

u  ( i ) 

- 

u  (i-1 

due 

= 

u(i-l) 

- 

u  ( i-2 

dud 

= 

dub 

- 

due 

due 

= 

dua 

- 

dub 

dp  a 

= 

P  (i  +  1) 

- 

P  ( i ) 

dpb 

= 

P  ( i ) 

- 

P  ( i-1 

dpc 

= 

P  ( i  —  1 ) 

- 

P  (i-2 

dpd 

= 

dpb 

- 

dpc 

dpe 

dp  a 

- 

dpb 

dza 

= 

z  (i  +  1) 

- 

z  (i) 

dzb 

= 

z  (i) 

- 

z  (i-1 

dze 

= 

z  ( i  —  1 ) 

- 

z  ( i-2 

dzd 

= 

dzb 

- 

dze 

dze 

= 

dza 

- 

dzb 

c  Check  monotonicity 


imon  =  1 

if 

(dra*drb 

.It. 

OdO) 

imon 

if 

(drb*drc 

•  It. 

OdO) 

imon 

if 

(dua*dub 

•  It. 

OdO) 

imon 

if 

(dub* due 

•  It. 

OdO) 

imon 

if 

(dpa*dpb 

•  It. 

OdO) 

imon 

if 

(dpb* dpc 

•  It. 

OdO) 

imon 

if 

(dza*dzb 

•  It. 

OdO) 

imon 

if 

(dzb*dzc 

•  It. 

OdO) 

imon 

if 

(imon  .eq. 

.  0) 

then 
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c  First-order  interface  is  non-monotonic 


rl 

=  r  ( i  —  1 

ul 

=  u  ( i  —  1 

Pi 

=  P ( i-1 

zl 

=  z  (i-1 

rr 

=  r  (i) 

ur 

=  u(i) 

pr 

=  P  ( i ) 

zr 

=  z  (i) 

else 

c  First-order 


& 


interface  is  monotonic 

denm  =  drb*drb  +  drc*drc  +  eps 
phi  =  (drb*drd  +  eps) /denm 
vhi  =  (drc*drd  +  eps) /denm 
rl  =  r(i-l)  +  iext* (epsm*phi*drc 

+  epsp*vhi*drb) 


& 


denm 

phi 

vhi 

rr 


dra*dra  +  drb*drb  +  eps 
(drb*dre  +  eps) /denm 
(dra*dre  +  eps) /denm 
r(i)  -  iext* (epsm*phi*dra 
+  epsp*vhi*drb) 


& 


denm 

phi 

vhi 

ul 


dub*dub  +  duc*duc  +  eps 
(dub*dud  +  eps) /denm 
(duc*dud  +  eps) /denm 
u(i-l)  +  iext* (epsm*phi*duc 
+  epsp*vhi*dub) 


& 


& 


& 


& 


denm 

phi 

vhi 

ur 


dua*dua  +  dub*dub  +  eps 
(dub*due  +  eps) /denm 
(dua*due  +  eps) /denm 
u(i)  -  iext* (epsm*phi*dua 
+  epsp*vhi*dub) 


denm 

phi 

vhi 

pi 


dpb*dpb  +  dpc*dpc  +  eps 
(dpb*dpd  +  eps) /denm 
(dpc*dpd  +  eps) /denm 
p(i-l)  +  iext* (epsm*phi*dpc 
+  epsp*vhi*dpb) 


denm 

phi 

vhi 

pr 


dpa*dpa  t  dpb*dpb  +  eps 
(dpb*dpe  +  eps) /denm 
(dpa*dpe  +  eps) /denm 
p(i)  -  iext* (epsm*phi*dpa 
+  epsp*vhi*dpb) 


denm 

phi 

vhi 

zl 


dzb*dzb  +  dzc*dzc  +  eps 
(dzb*dzd  +  eps) /denm 
(dzc*dzd  +  eps) /denm 
z(i-l)  +  iext* (epsm*phi*dzc 
+  epsp*vhi*dzb) 
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denm  = 

dza*dza 

+  dzb* 

dzb  +  eps 

phi  = 

(dzb*dze 

+  eps 

)  /denm 

vhi  = 

(dza*dze 

+  eps 

)  /denm 

zr  = 

z  (i) 

iext* 

(epsm*phi*dza 

&  +  epsp*vhi*dzb) 

endif 


else 

write ( * , * ) 
write  (  * , * ) 
write  (  * , * ) 
endif 


I  I 

'  Unknown  limiting  strategy' 

I  I 


c  Set  ceiling  on  zl,  zr 


zl 

= 

min 

zr 

= 

min 

else 

c  First 

order  at 

imax 

rl  = 

r 

(i-l 

ul  = 

u 

(i-1 

pi  = 

P 

(i-l 

zl  = 

z 

(i-1 

rr  = 

r 

(i) 

ur  = 

u 

(i) 

pr  = 

P 

(i) 

zr  = 

z 

(i) 

endif 

c 

zzl  (i) 

= 

zl 

c 

zzr (i) 

= 

zr 

ldO ) 
ldO ) 


c  Final  monotonicity  check 


imon  = 

0 

rat 

= 

(r  (i) 

-  r (i-1) ) * (rr 

-  rl ) 

if 

(rat 

.It. 

OdO)  imon  =  1 

rat 

= 

(u(i) 

-  u (i-1) ) * (ur 

-  ul) 

if 

(rat 

.It. 

OdO)  imon  =  2 

rat 

= 

(p  ( i ) 

-  p (i-1) ) * (pr 

-  pl) 

if 

(rat 

.It. 

OdO)  imon  =  3 

rat 

= 

(z  (i) 

-  z (i-1) ) * (zr 

-  zl) 

if 

(rat 

.It. 

OdO)  imon  =  4 

c  Set  first  order  interface 

if  (imon  .ne.  0)  then 
rl  =  r(i-l) 
rr  =  r ( i ) 
ul  =  u ( i — 1 ) 
ur  =  u  ( i ) 
pi  =  p(i-l) 
pr  =  p (i) 
zl  =  z  (i-1) 
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zr 


z  (i) 


endif 

c  Calculate  internal  energy 

if  (ieos  .eq.  0)  then 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  CPG  EOS 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
el  =  pl/gaml/rl  -  zl*qdet0 

zr*qdet0 


er 


=  pr/gaml/rr  - 
else  if  (ieos  .eq.  1)  then 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  JWL  EOS 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


rht 

=  rl/rO 

rhti 

=  ldO/rht 

ri 

=  IdO/rl 

tmp 

=  pi  -  aj*(ld0  -  wrl*rl) 
-  bj*(ld0  -  wr2*rl) 

*dexp 

*dexp 

el 

=  tmp*ri/wj  -  zl*qdet0 

rht 

=  rr/rO 

rhti 

=  ldO/rht 

ri 

=  IdO/rr 

tmp 

=  pr  -  aj*(ld0  -  wrl*rr) 
-  bj*(ld0  -  wr2*rr) 

*dexp 

*dexp 

er 

=  tmp*ri/wj  -  zr*qdet0 

-rh2*ri) 


-rh2*ri) 


else  if  (ieos  .eq.  2)  then 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Hayes-I/ JWL  EOS 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Left  of  interface;  set  arguments 
ra  =  rl 
za  =  zl 
rz  =  ra*za 
omz  =  ldO  -  za 

c  Solid  phase  limit 

if  (za  .le.  ztoll)  then 

el  =  pl/gh  +  beta*r0/ra  +  t4* ( (ra/rO ) **alf a)  -  t7 

c  Mixed  phases 

else  if  (ztoll  .It.  za  .and.  za  .It.  ztol2)  then 

c  Evaluate  denominator  function 

bot  =  omz/gh  +  ldO/wj/ra 
if  (bot  .It.  Id-10)  then 
write  (  * , * )  '  ' 

write  (*,*)  '  Zero  denonimator  term.' 
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write ( * , * )  '  ' 

stop 
endif 

c  Evaluate  numerator  functions 

top(2)  =  omz  -  rO/ra 

top(3)  =  (omz**nh) * ( (ra/rO ) **alfa) 

top (4)  =  omz 

top(5)  =  (ldO/wj/ra  -  za/rhl ) *dexp ( -rhl /rz ) 

top(6)  =  (ldO/wj/ra  -  za/rh2 ) *dexp ( -rh2 /rz ) 

top(7)  =  za 

el  =  bot*pl 
do  nn  =  2,7 

el  =  el  -  c (nn) *top (nn) 
enddo 

c  Gas  phase  limit 
else 

el  =  pl/wj/ra 

-  a j * ( IdO /wj /ra  -  ldO/rhl ) *dexp ( -rhl /ra) 

-  b j * ( IdO /wj /ra  -  IdO /rh2 ) *dexp ( -rh2 /ra) 

-  qdetO  -  eO 

endif 

c  Right  of  interface;  set  arguments 
ra  =  rr 
za  =  zr 
rz  =  ra*za 
omz  =  IdO  -  za 

c  Solid  phase  limit 

if  (za  .le.  ztoll)  then 

er  =  pr/gh  +  beta*rO/ra  +  t4* ( (ra/rO ) **alf a)  -  t7 

c  Mixed  phases 

else  if  (ztoll  .It.  za  .and.  za  .It.  ztol2)  then 

c  Evaluate  denominator  function 

bot  =  omz/gh  +  ldO/wj/ra 
if  (bot  .It.  Id-10)  then 
write  (  * , * )  '  ' 

write (*,*)  '  Zero  denonimator  term.' 

write  (  * , * )  '  ' 

stop 
endif 

c  Evaluate  numerator  functions 

top(2)  =  omz  -  rO/ra 

top(3)  =  (omz**nh) *( (ra/rO ) **alfa) 

top (4)  =  omz 

top(5)  =  (ldO/wj/ra  -  za/rhl ) *dexp ( -rhl /rz ) 

top(6)  =  (ldO/wj/ra  -  za/rh2 ) *dexp ( -rh2/rz ) 

top(7)  =  za 
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& 

& 
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c  Compute  internal  energy 
er  =  bot*pr 
do  nn  =  2,7 

er  =  er  -  c (nn) *top (nn) 
enddo 

c  Gas  phase  limit 
else 

er  =  pr/wj/ra 

-  a j * ( ldO/wj /ra  -  ldO /rhl ) *dexp ( -rhl /ra) 

-  b j * ( ldO/wj /ra  -  Id0/rh2 ) *dexp ( -rh2 /ra) 

-  qdetO  -  eO 

endif 

else  if  (ieos  .eq.  3)  then 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Hayes-II/JWL  EOS 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Left  of  interface;  set  arguments 
ra  =  rl 
za  =  zl 
rz  =  ra*za 
omz  =  ldO  -  za 

c  Solid  phase  limit 

if  (za  .le.  ztoll)  then 

el  =  pl/gh  +  beta*rO/ra  +  t4* ( (ra/rO ) **alf a)  -  t7 
&  -  hi /gh/nh* (( (ra/rO ) **nh)  -  ldO) 

c  Mixed  phases 

else  if  (ztoll  .It.  za  .and.  za  .It.  ztol2)  then 

c  Evaluate  denominator  function 

bot  =  omz/gh  +  ldO/wj/ra 
if  (bot  .It.  Id-10)  then 
write  (  * , * )  '  ' 

write  (*,*)  '  Zero  denonimator  term.' 

write  (  * , * )  '  ' 

stop 
endif 

c  Evaluate  numerator  functions 

top(2)  =  omz  -  rO/ra 

top(3)  =  (omz**nh) *( (ra/rO ) **alfa) 

top (4)  =  omz 

top(5)  =  (ldO/wj/ra  -  za/rhl ) *dexp ( -rhl /rz ) 

top(6)  =  (ldO/wj/ra  -  za/rh2 ) *dexp ( -rh2/rz ) 

top(7)  =  za 

top(8)  =  (omz**nhpl )*( (ra/rO ) **nh)  +  za  -  ldO 

el  =  bot*pl 
do  nn  =  2,8 
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& 

& 

& 


el  =  el  -  c (nn) *top (nn) 
enddo 

c  Gas  phase  limit 
else 

el  =  pl/wj/ra 

&  -  a j * ( ldO /wj /ra  -  ldO /rhl ) *dexp ( -rhl /ra) 

&  -  b j * ( ldO /wj /ra  -  ldO /rh2 ) *dexp ( -rh2/ra) 

&  -  qdetO  -  eO 

endif 

c  Right  of  interface;  set  arguments 


ra 

=  rr 

za 

=  zr 

rz 

=  ra* 

omz 

=  ldO 

c  Solid  phase  limit 

if  (za  .le.  ztoll)  then 

er  =  pr/gh  +  beta*rO/ra  +  t4* ( (ra/rO ) **alf a)  -  t7 
&  -  hi /gh/nh* ((( ra/rO ) **nh)  -  ldO) 

c  Mixed  phases 

else  if  (ztoll  .It.  za  .and.  za  .It.  ztol2)  then 

c  Evaluate  denominator  function 

bot  =  omz/gh  +  ldO/wj/ra 
if  (bot  .It.  Id-10)  then 
write  (  * , * )  '  ' 

write (*,*)  '  Zero  denonimator  term.' 

write  (  * , * )  '  ' 

stop 
endif 

c  Evaluate  numerator  functions 

top(2)  =  omz  -  rO/ra 

top(3)  =  (omz**nh) *( (ra/rO ) **alfa) 

top (4)  =  omz 

top(5)  =  (ldO/wj/ra  -  za/rhl ) *dexp ( -rhl /rz ) 

top(6)  =  (ldO/wj/ra  -  za/rh2 ) *dexp ( -rh2/rz ) 

top(7)  =  za 

top(8)  =  (omz**nhpl )*( (ra/rO ) **nh)  +  za  -  ldO 

c  Compute  internal  energy 
er  =  bot*pr 
do  nn  =  2,8 

er  =  er  -  c (nn) *top (nn) 
enddo 

c  Gas  phase  limit 
else 

er  =  pr/wj/ra 

&  -  a j * ( ldO /wj /ra  -  ldO /rhl ) *dexp ( -rhl /ra) 
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&  -  b j * ( IdO /wj /ra  -  ldO /rh2 ) *dexp ( -rh2/ra) 

&  -  qdetO  -  eO 

endif 

else 

write  ( * , * )  '  ' 

write  (*,*)  '  Unknown  EOS' 

write  (  * , * )  '  ' 

stop 
endif 


c  Total  energy/mass 

eel  =  el  +  0.5d0*ul*ul 
hhl  =  eel  +  pl/rl 

eer  =  er  +  0.5d0*ur*ur 
hhr  =  eer  +  pr/rr 


c 

if  (imon  .ne. 

0)  then 

c 

write ( * , * )  ' 

1 

c 

write  (  * , * )  ' 

Monotonicity  violation 

c 

write ( * , * )  ' 

i  = 

i 

c 

write ( * , * )  ' 

I 

c 

write ( * , * )  ' 

r  ( i  —  1 ) 

=  '  ,  r  ( i  —  1 ) 

c 

write ( * , * )  ' 

rl 

=  '  ,  rl 

c 

write  (  * , * )  ' 

rr 

=  '  ,  rr 

c 

write ( * , * )  ' 

r  (i) 

=  '  ,  r (i) 

c 

write ( * , * )  ' 

i 

c 

write ( * , * )  ' 

u(i-l) 

=  ' , u ( i— 1 ) 

c 

write  (  * , * )  ' 

ul 

i — 1 

3 

II 

c 

write ( * , * )  ' 

ur 

=  '  ,  ur 

c 

write ( * , * )  ' 

u(i) 

=  '  ,  u  (i) 

c 

write ( * , * )  ' 

I 

c 

write ( * , * )  ' 

P(i-l) 

=  '  ,P ( i  —  1 ) 

c 

write  (  * , * )  ' 

pi 

=  ',pl 

c 

write ( * , * )  ' 

pr 

=  '  ,pr 

c 

write ( * , * )  ' 

P  ( i ) 

=  ',P(D 

c 

write ( * , * )  ' 

1 

c 

pause 

c 

endif 

c80 

format (2x,dl2.6,2x 

,  dl2 . 6, 

2x, dl2 . 6) 

c 

if  (n  .eq.  177 

)  then 

c 

write (25, 80) 

r  ( i  —  1 ) 

,  rr, r (i) 

c 

endif 

c  Roe 

averages 

if  (iav  .eq.  1 

)  then 

sqrl  =  dsqrt(rl) 
sqrr  =  dsqrt(rr) 
rsumi  =  ldO/ (sqrl  +  sqrr) 


'  ,  imon 


rav 

=  sqrl*sqrr 

uav 

=  (sqrl*ul  + 

sqrr*ur) *rsumi 

zav 

=  (sqrl*zl  + 

sqrr*zr) *rsumi 

zav 

=  min ( zav,  ldO 

) 
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(sqrl*el  +  sqrr*er) *rsumi 
(sqrl*hhl  +  sqrr*hhr) *rsumi 


eav 
hav 
else 

c  Test  arithmetic  averages 

rav  =  0.5d0*(rl  +  rr) 

uav  =  0.5d0*(ul  +  ur) 

zav  =  0.5d0*(zl  +  zr) 

eav  =  0.5d0*(el  +  er) 

hav  =  0.5d0*(hhl  +  hhr) 
endif 

pav  =  rav* (hav  -  eav  -  0 . 5d0*uav*uav) 

c  Calculate  averaged  pressure  derivatives 
if  (ieos  .eq.  0)  then 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c  CPG  EOS 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
dpdr  =  gaml*eav  +  gaml*zav*qdetO 
dpde  =  gaml*rav 

dpdz  =  gaml*rav*qdetO 

else  if  (ieos  .eq.  1)  then 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c  JWL  EOS 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ri  =  IdO/rav 

tmp  =  a j * (rhl*ri*ri  -  wj*ri  -  wj /rhl ) *dexp ( -rhl*ri ) 

tmp  =  tmp  +  bj * (rh2*ri*ri  -  wj*ri  -  wj /rh2 ) *dexp ( -rh2*ri ) 

dpdr  =  tmp  +  wj*eav  +  wj*zav*qdetO 

dpde  =  wj*rav 

dpdz  =  wj*rav*qdetO 

else  if  (ieos  .eq.  2)  then 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Hayes-I/ JWL  EOS 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


ra 

=  rav 

ra2 

=  ra*: 

ea 

=  eav 

za 

=  zav 

rz 

=  ra* 

omz 

=  ldO 

c  Solid  phase  limit 

if  (za  .le.  ztoll)  then 

dpdr  =  beta*r0 *gh/ra2  -  alfa*gh*t4* (ra** (alfa-ldO)  ) 
&  /  r0**alfa 

dpdz  =  gh*ea  -  beta*rO*gh/ra  +  alfa*gh*t4 

&  *  ( (ra/rO ) **alf a) 
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dpde  =  gh 


c  Mixed  phases 

else  if  (ztoll  .It.  za  .and.  za  .It.  ztol2)  then 

c  Evaluate  denominator  functions 

bot  =  omz/gh  +  ldO/wj/ra 
if  (bot  .It.  Id-10)  then 
write  (  * , * )  '  ' 

write (*,*)  '  Zero  denonimator  term.' 

write  (  * , * )  '  ' 

stop 
endif 

bot2  =  bot*bot 
botr  =  -Id0/wj/ra2 
botz  =  -ldO/gh 

c  Evaluate  numerator  functions 
top(l)  =  ea 
top(2)  =  omz  -  rO/ra 
top(3)  =  (omz**nh) * ( (ra/rO ) **alfa) 
top (4)  =  omz 

top(5)  =  (ldO/wj/ra  -  za/rhl ) *dexp ( -rhl /rz ) 
top(6)  =  (ldO/wj/ra  -  za/rh2 ) *dexp ( -rh2 /rz ) 
top(7)  =  za 


c  Compute  derivatives  for  numerator  functions 


topr ( 1 ) 

=  OdO 

topr (2 ) 

=  r0/ra2 

topr ( 3 ) 

=  alfa/rO* (omz**nh) * ( (ra/rO) ** (alfa-ldO) ) 

topr ( 4 ) 

=  OdO 

topr ( 5 ) 

=  (rhl/wj/rz  -  ldO/wj  -  ldO ) *dexp ( -rhl /rz ) /ra2 

topr ( 6 ) 

=  (rh2/wj/rz  -  ldO/wj  -  ldO ) *dexp ( -rh2 /rz ) /ra2 

topr ( 7 ) 

=  OdO 

topz ( 1 ) 

=  OdO 

topz (2 ) 

=  -ldO 

topz (3) 

=  -nh* ( (omz*ra/r0 ) **alf a) 

topz ( 4 ) 

=  -ldO 

topz  ( 5 ) 

=  (rhl/wj /rz/rz  -  ldO/rz  -  ldO/rhl) *dexp (- 

rhl / rz ) 

topz  ( 6 ) 

=  (rh2/wj /rz/rz  -  ldO/rz  -  ldO / rh2 ) *dexp ( - 

rh2/ rz ) 

topz ( 7 ) 

=  ldO 

c  Compute  density  and  internal  energy  derivatives  of  pressure 

dpdr  = 

OdO 

dpdz  = 

OdO 

do  nn  = 

=  1,7 

dpdr 

=  dpdr  +  c (nn) * (bot*topr (nn)  -  botr*top (nn) ) 

dpdz 

=  dpdz  +  c (nn) * (bot*topz  (nn)  -  botz*top  (nn) ) 

enddo 

dpdr  = 

dpdr/bot2 

dpdz  = 

dpdz/bot2 

dpde  = 

ldO /bot 
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c  Gas  phase  limit 
else 


& 

& 

& 


& 

& 


dpdr  =  wj  *ei ( i ) 

+  aj*(rhl/ra2  -  w j /ra  -  wj /rhl ) *dexp ( -rhl /ra) 
+  bj*(rh2/ra2  -  w j /ra  -  wj /rh2 ) *dexp ( -rh2 /ra) 
+  w j  * (qdetO  +  eO ) 

dpdz  =  aj*(rhl/ra  -  wj  -  wj *ra/rhl ) *dexp ( -rhl /ra) 

+  bj*(rh2/ra  -  wj  -  wj *ra/rh2 ) *dexp ( -rh2 /ra) 

+  ra*wj* (qdetO  +  eO) 

dpde  =  wj*ra 

endif 


else  if  (ieos  .eq.  3)  then 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Hayes-II/JWL  EOS 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


ra 

=  rav 

ra2 

=  ra*. 

ea 

=  eav 

za 

=  zav 

rz 

=  ra* 

omz 

=  ldO 

c  Solid  phase  limit 

if  (za  .le.  ztoll)  then 


& 

& 


dpdr  =  beta*rO *gh/ra2  -  alfa*gh*t4* (ra** (alfa-ldO) ) 
/  rO**alfa 

+  hl/rO* ( (ra/rO ) **nhml ) 


& 

& 


dpdz  =  gh*ea  -  beta*rO*gh/ra  +  alfa*gh*t4 
*  ( (ra/rO ) **alf a) 

+  hl/nh*(ldO  -  nhpl* ( (ra/rO ) **nh) ) 
dpde  =  gh 


c  Mixed  phases 

else  if  (ztoll  .It.  za  .and.  za  .It.  ztol2)  then 


c  Evaluate  denominator  functions 

bot  =  omz/gh  +  ldO/wj/ra 
if  (bot  .It.  Id-10)  then 
write  (  * , * )  '  ' 

write (*,*)  '  Zero  denonimator  term.' 

write  (  * , * )  '  ' 

stop 
endif 

bot2  =  bot*bot 
botr  =  -Id0/wj/ra2 
botz  =  -IdO/gh 

c  Evaluate  numerator  functions 
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top  (1) 

=  ea 

top (2) 

=  omz  -  rO/ra 

top  ( 3 ) 

=  (omz**nh) *( (ra/rO ) ** 

alfa) 

top (4) 

=  omz 

top ( 5 ) 

=  (IdO/wj/ra  -  za/rhl) 

*dexp ( -rhl /rz ) 

top ( 6 ) 

=  (IdO/wj/ra  -  za/rh2) 

*dexp ( -rh2 /rz ) 

top (7) 

=  za 

top ( 8 ) 

=  (omz**nhpl )*( (ra/rO  ) 

**nh)  +  za  -  ldO 

c  Compute  derivatives  for  numerator  functions 


rhl / rz ) 
rh2/ rz) 


topr 

(1) 

= 

OdO 

topr 

(2) 

= 

rO / ra2 

topr 

(3) 

= 

alfa/ rO* (omz 

topr 

(4) 

= 

OdO 

topr 

(5) 

= 

(rhl/wj/rz  - 

topr 

(6) 

= 

(rh2/wj/rz  - 

topr 

(7) 

= 

OdO 

topr 

(8) 

= 

nh/ rO* (omz**i 

topz 

(1) 

= 

OdO 

topz 

(2) 

= 

-ldO 

topz 

(3) 

= 

-nh* ( (omz*ra 

topz 

(4) 

= 

-ldO 

topz 

(5) 

= 

(rhl/wj /rz/r 

topz 

(6) 

= 

(rh2/wj /rz/r 

topz 

(7) 

= 

ldO 

topz 

(8) 

= 

ldO  -  nhpl* ( 

iity 

and 

.  internal  energy 

dpdr 

= 

OdO 

dpdz 

= 

OdO 

do  nn  =  1,8 

dpdr  =  dpdr  +  c (nn) * (bot*topr (nn) 
dpdz  =  dpdz  +  c  (nn) * (bot*topz (nn) 
enddo 

dpdr  =  dpdr/bot2 
dpdz  =  dpdz/bot2 
dpde  =  ldO/bot 


botr*top (nn) ) 
botz*top (nn) ) 


c  Gas  phase  limit 
else 


dpdr  =  wj *ei ( i ) 

+  aj*(rhl/ra2  -  w j /ra  -  wj /rhl ) *dexp ( -rhl /ra) 
+  bj*(rh2/ra2  -  w j /ra  -  wj / rh2 ) *dexp ( -rh2 /ra) 
+  w j  * (qdetO  +  eO ) 

dpdz  =  aj*(rhl/ra  -  wj  -  wj *ra/rhl ) *dexp ( -rhl /ra) 

+  bj*(rh2/ra  -  wj  -  wj *ra/rh2 ) *dexp ( -rh2 /ra) 

+  ra*wj* (qdetO  +  eO) 

dpde  =  wj*ra 

endif 
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else 

write  ( * , * 

)  '  ' 

write  ( * , * 

)  '  Unknown  EOS ' 

write  ( * , * 

)  '  ' 

stop 

endif 

c  Calculate  averaged 

speed  of 

sound 

if  (dpdr  .It.  OdO) 

dpdr  =  dabs (dpdr) 

a2  =  dpdr 

+  pav*dpde/rav/rav 

if  (a2  .It. 

OdO)  then 

write ( * , * 

)  '  a2  < 

0  !  ' 

write ( *  ,  * 

)  '  i  = 

M 

write ( * , * 

)  '  eav 

=  ' , eav 

write ( * ,  * 

)  '  el 

=  \el 

write  ( * , * 

)  '  er 

=  '  ,  er 

write  ( * , * 

)  '  rav 

=  ' , rav 

write  ( * , * 

)  '  pav 

=  '  ,  pav 

write  ( * , * 

)  '  pl 

=  -,pl 

write  ( * , * 

)  '  pr 

=  '  ,pr 

write  ( * , * 

)  '  zav 

=  ' , zav 

write  ( * , * 

)  '  dpdr 

=  ' , dpdr 

write  ( * , * 

)  '  dpde 

=  ' , dpde 

write  ( * , * 

)  '  ' 

write  ( * , * 

)  '  r+1 

=  ' , rp ( i ) 

write  ( * , * 

)  '  u+1 

=  ' r up ( i ) 

write  ( * , * 

)  '  P+1 

=  ' f  PP ( i ) 

write  ( * , * 

)  '  z  +  1 

=  '  ,  zp (i) 

write  ( * , * 

)  '  ' 

write  ( * , * 

t — 1 

1 

u 

=  '  ,  rp (i-1) 

write  (  * , * 

t — 1 

1 

3 

=  ' , up ( i-1 ) 

write  ( * , * 

)  '  p-i 

=  ' rPP (i-1) 

write  ( * , * 

)  '  z-l 

=  ' , zp (i-1) 

write  ( * , * 

)  '  ' 

stop 

endif 

aav  =  dsqrt (a2 ) 

if  (idbgf  . 

eg.  1)  then 

write  (  * , * ) 

'  rav  = 

'  ,  rav 

write  (  * , * ) 

'  uav  = 

'  ,  uav 

write ( * , * ) 

’  zav  = 

'  ,  zav 

write ( * , * ) 

'  eav  = 

'  ,  eav 

write  (  * , * ) 

'  hav  = 

’  ,  hav 

write  (  * , * ) 

'  pav  = 

'  r  PaV 

write ( * , * ) 

'  aav  = 

'  ,  aav 

write ( * , * ) 

?  i 

endif 

c  Eigenvalues 

aeg(l)  =  dabs (uav  - 

aav) 

aeg(2)  =  dabs (uav) 

aeg(3)  =  dabs (uav) 

aeg(4)  =  dabs (uav  + 

aav) 
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if  (ldbgf  .eq.  1)  then 


write ( * , * ) 

'  aegl 

write ( * , * ) 

'  aeg2 

write  (  * , * ) 

'  aeg3 

write  (  * , * ) 

'  aeg4 

write ( * , * ) 

I  1 

pause 

endif 

c  Right  eigenvectors 


evr  ( 1 , 

1) 

= 

ldO 

evr  ( 1 , 

2) 

= 

ldO 

evr  ( 1 , 

3) 

= 

ldO 

evr  ( 1 , 

4) 

= 

ldO 

evr (2 , 

1) 

= 

uav 

- 

aav 

evr (2 , 

2) 

= 

uav 

evr (2 , 

3) 

= 

uav 

evr (2 , 

4) 

= 

uav 

+ 

aav 

evr ( 3 , 

1) 

= 

hav 

_ 

uav 

evr ( 3 , 

2) 

= 

hav 

- 

rav 

evr ( 3 , 

3) 

= 

hav 

- 

rav 

evr ( 3 , 

4) 

= 

hav 

+ 

uav 

evr ( 4 , 

1) 

= 

zav 

evr  ( 4 , 

2) 

= 

OdO 

evr  ( 4 , 

3) 

= 

ldO 

evr  ( 4 , 

4) 

= 

zav 

, aeg  (1) 
, aeg  (2) 
, aeg (3) 
, aeg  (4) 


aav 

a2/dpde  +  zav*dpdz/dpde 
a2/dpde  +  (zav  -  ldO ) *dpdz/dpde 
aav 


if  (idbgf  .eq.  1)  then 
write  (*, *)  ' EVR:  ' 

write  (*,71)  evr (1,1), evr (1,2) , evr (1,3) , evr (1,4) 
write ( * , 71)  evr (2,1), evr (2,2) , evr (2,3) , evr (2,4) 
write (*,71)  evr (3,1), evr (3,2) , evr (3,3) , evr (3,4) 
write (*,71)  evr (4,1), evr (4,2) , evr (4,3) , evr (4,4) 
write ( * , * )  '  ' 

endif 


c  I  R  I 

detr  =  -2d0*rav*a2*aav/dpde 

c  Compute  primitive  variables  differences 
delr  =  rr  -  rl 
delv  =  ur  -  ul 
delp  =  pr  -  pi 
delz  =  zr  -  zl 


c  Compute  characteristic  wave  magnitudes 
omz  =  ldO  -  zav 

cwm ( 1 )  =  cl2* (delp/aav/aav  -  rav*delv/aav) 
cwm(2)  =  omz* (delr  -  delp/aav/aav)  -  rav*delz 
cwm(3)  =  zav* (delr  -  delp/aav/aav)  +  rav*delz 
cwm(4)  =  cl2* (delp/aav/aav  +  rav*delv/aav) 


c  Compute  R  eg  L  dq 
do  1  =  1,4 
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vn ( 1 )  =  OdO 
do  m  =  1,4 

vn(l)  =  vn(l)  +  aeg(m) *cwm(m) *evr(l,m) 
enddo 
enddo 


if  (idbgf 
write ( * , * ) 
write  (  * , * ) 
write  (  * , * ) 
write ( * , * ) 
write ( * , * ) 
endif 


. eq.  1)  then 
'  vnl  =  '  ,  vn ( 1 ) 
'  vn2  =  ' , vn ( 2 ) 
'  vn3  =  '  ,  vn ( 3 ) 
'  vn4  =  ' , vn ( 4 ) 


c  Compute  the  Euler  flux 


fl  (1) 

=  rl*ul 

fl  (2) 

=  rl*ul*ul  + 

Pi 

fl  (3) 

=  rl*ul*hhl 

fl  (4) 

=  rl*ul*zl 

fr  (1) 

=  rr*ur 

fr  (2) 

=  rr*ur*ur  + 

pr 

fr  (3) 

=  rr*ur*hhr 

f  r  ( 4  ) 

=  rr*ur*zr 

c  Compute  numerical  flux 
do  1  =  1,4 

fn(i,l)  =  0.5d0*(fl(l)  +  fr(l)  -  vn(l)) 
enddo 

if  (idbgf  .eg.  1)  then 


write ( * , * ) 

1 

FL:  ' 

write ( * , * ) 

1 

fll  = 

'  ,  fl (1) 

write  (  * , * ) 

1 

f  12  = 

'  ,  fl (2) 

write  (  * , * ) 

1 

f  13  = 

'  ,  fl (3) 

write ( * , * ) 

1 

f  14  = 

'  ,  fl (4 ) 

write ( * , * ) 

1 

1 

write  (  * , * ) 

1 

FR:  ' 

write  (  * , * ) 

1 

f  rl  = 

'  rfr (1) 

write ( * , * ) 

1 

f  r2  = 

'  /  fr  (2) 

write ( * , * ) 

1 

f  r3  = 

' ;  fr (3) 

write ( * , * ) 

1 

f  r4  = 

' / fr  (4) 

write  (  * , * ) 

1 

1 

write ( * , * ) 

I 

FN:  ' 

write ( * , * ) 

1 

fnl  = 

' , fn  (i,  1 

write ( * , * ) 

1 

fn2  = 

',fn(i,2 

write  (  * , * ) 

1 

f  n3  = 

'  ,  fn (i, 3 

write  (  * , * ) 

1 

fn4  = 

',fn(i,4 

write  ( * , * ) 
pause 
endif 
enddo 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  End  of  flux  calculation  loop 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Advance  the  solution  in  time 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

do  i  =  1 , imax-1 
do  1  =  1,4 

dqv(l)  =  dt/dx* (fn (i  +  1 , 1)  -  fn(i,l)) 

enddo 

do  1  =  1,4 

qvp(i,l)  =  qv(i,l)  -  dqv(l)  +  dt*s(i,l) 
enddo 
enddo 

c  Extract  primitive  variables 
do  i  =  1 , imax-1 


rp  ( i ) 

qvp ( i , 

1) 

up  ( i ) 

qvp ( i , 

2)  / qvp  ( i 

,D 

etp(i)  = 

qvp ( i , 

3) /qvp (i 

,D 

zp  ( i ) 

qvp ( i , 

4) / qvp  ( i 

,D 

zp  ( i ) 

min  ( zp 

>(i) , ldO) 

zp  ( i ) 

max ( zp 

> ( i ) , OdO ) 

if  ( zp ( i ) 

1  .  It . 

ld-99) 

zp  (i) 

=  OdO 

if  (zp(i; 

i  .  ge . 

0 . 99d0 ) 

zp  (i) 

=  ldO 

eip(i)  = 

etp (i) 

-  0 . 5d0 

*up (i) 

*up  ( i 

tk(i) 

tkO 

if  (rp (i; 

i  .  le . 

OdO)  then 

write  (  * , * )  ' 

? 

write(*,*)  ' 

Negative/Zero  density' 

write (*,*)  ' 

i  =  '  ,  i 

write(*,*)  ' 

r  =  ' , rp ( i ) 

write(*,*)  ' 

u  =  ' , up (i) 

write ( * , * )  ' 

e  =  ' , etp (i) 

write(*,*)  ' 

z  =  ' , zp (i) 

write(*,*)  ' 

1 

write(*,*)  ' 

Program  STOP' 

write (*,*)  ' 

stop 
endif 

1 

c  If 

internal  energy  is 

negative,  apply  a  fix 

if  (eip(i)  . le 

.  OdO)  then 

c 

write(*,*)  ' 

1 

c 

write (*,*)  ' 

Negative/Zero  internal  energy 

c 

write(*,*)  ' 

i  =  '  ,  i 

c 

write(*,*)  ' 

r  =  ' , rp ( i ) 

c 

write  (  * , * )  ' 

u  =  ' , up ( i ) 

c 

write(*,*)  ' 

E  =  '  , etp ( i ) 

c 

write (*,*)  ' 

e  =  ' , eip ( i ) 

c 

write(*,*)  ' 

z  =  '  ,  zp (i) 
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c 

c 


write  (  * , * )  '  ' 

write  ( * , * )  '  ief x 


'  ,  iefx 


if  (iefx  .eq.  0)  then 
c  Absolute  value  |e|  fix 

eip(i)  =  dabs(eip(i)) 

else  if  (iefx  .eq.  1)  then 
c  Pressure  estimation  fix 
c  Estimate  pressure  using  JWL  EOS 

pest  =  a j *dexp ( -rhl /rp ( i ) )  +  bj *dexp (-rh2/rp (i) ) 

&  +  cjh* (rp (i) ** (ldO  +  wj ) ) 

c  Compute  detonation  e  based  on  JWL  pressure 
eip(i)  =  ldO /wj /rp ( i ) * 

&  (  cjh* (rp (i) ** (ldO  +  wj )  ) 

&  +  aj *wj *rp (i) /rhl*dexp (-rhl/rp (i) ) 

&  +  bj *wj *rp (i) /rh2*dexp (-rh2/rp (i) )  ) 

write ( * , * )  '  ' 

write (*,*)  '  pest  =  ' ,pest 
write(*,*)  '  eest  =  ',eip(i) 
pause 

else  if  (iefx  .eq.  2)  then 
c  Time-lagged  velocity  fix 

eip(i)  =  etp(i)  -  0 . 5d0 *u ( i ) *u ( i ) 

else  if  (iefx  .eq.  3)  then 
c  Zero  kinetic  energy  fix 

eip(i)  =  etp(i) 

else 

write  ( * , * )  '  ' 

write (*,*)  '  Unknown  iefx  value.' 

write  (  * , * )  '  ' 

stop 
endif 

c  pause 

endif 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Calculate  pressure  and  its  derivatives 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
if  (ieos  .eq.  0)  then 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  CPG  EOS 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
pp(i)  =  gaml*rp (i) *eip (i)  +  gaml*rp (i) *zp (i) *qdet0 

dpdr  =  gaml*eip(i)  +  gaml*zp (i) *qdet0 

dpde  =  garni *rp(i) 

dpdz  =  gaml*rp (i) *qdet0 

else  if  (ieos  .eq.  1)  then 
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c 

c 

c 

c 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  JWL  EOS 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
rht  =  rp ( i ) /rO 
rhti  =  ldO/rht 
ri  =  IdO/ rp  (i) 

tmp  =  aj*(ldO  -  wrl*rp (i) ) *dexp (~rhl*ri) 

tmp  =  tmp  +  bj*(ldO  -  wr2*rp (i) ) *dexp (-rh2*ri) 

pp(i)  =  tmp  +  wj *rp ( i ) *eip ( i )  +  wj *rp ( i ) * zp ( i ) *qdetO 

tmp  =  a j * (rhl*ri*ri  -  wj*ri  -  wj /rhl) *dexp (-rhl*ri) 

tmp  =  tmp  +  bj * (rh2*ri*ri  -  wj*ri  -  wj /rh2 ) *dexp ( -rh2*ri ) 

dpdr  =  tmp  +  wj*eip(i)  +  wj * zp ( i ) *qdetO 

dpde  =  wj*rp(i) 

dpdz  =  wj *rp ( i ) *qdetO 

else  if  (ieos  .eq.  2)  then 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Hayes-I/ JWL  EOS 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


ra 

=  rp ( i ) 

ra2 

=  ra*ra 

ea 

=  eip ( i 

za 

=  zp(i) 

rz 

=  ra*za 

omz 

=  IdO  - 

c  Solid  phase  limit 

if  (za  .le.  ztoll)  then 

c  Compute  pressure  and  its  derivatives 

pp(i)  =  gh* (ea  -  beta*rO/ra  -  t4* ( (ra/rO ) **alf a) 

&  +  t7 ) 

dpdr  =  beta*r0*gh/ra2  -  alfa*gh*t4* (ra** (alfa-ldO) ) 
&  / (rO**alfa) 

dpdz  =  gh*ea  -  beta*rO*gh/ra  +  alfa*gh*t4 

&  *  ( (ra/rO ) **alf a) 

dpde  =  gh 

c  Mixed  phases 

else  if  (ztoll  .It.  za  .and.  za  .It.  ztol2)  then 

c  Evaluate  denominator  functions 

bot  =  omz/gh  +  ldO/wj/ra 
if  (bot  .It.  Id-10)  then 
write  (  * , * )  '  ' 

write (*,*)  '  Zero  denonimator  term.' 

write  (  * , * )  '  ' 

stop 
endif 
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bot2 

botr 

botz 


=  bot*bot 
=  -Id0/wj/ra2 
=  -IdO/gh 


numerator 

functions 

top  (1) 

=  ea 

top  (2) 

=  omz  -  rO/ra 

top ( 3 ) 

=  (omz**nh) * ( (ra/rO ) ** 

alfa) 

top (4) 

=  omz 

top ( 5 ) 

=  (IdO/wj/ra  -  za/rhl) 

*dexp ( -rhl / rz 

top ( 6 ) 

=  (IdO/wj/ra  -  za/rh2) 

*dexp ( -rh2 / rz 

top  (7) 

=  za 

c  Compute  derivatives 

for  numerator  functions 

topr  ( 1 ) 

= 

OdO 

topr (2 ) 

= 

r0/ ra2 

topr  ( 3 ) 

= 

alfa/rO* (omz**nh) * ( (ra/rO) ** (alfa-ldO) ) 

topr ( 4 ) 

= 

OdO 

topr ( 5 ) 

= 

(rhl/wj/rz  -  ldO/wj  -  IdO ) *dexp ( -rhl /rz ) 

topr ( 6 ) 

= 

(rh2/wj/rz  -  ldO/wj  -  IdO ) *dexp ( -rh2/rz ) 

topr  ( 7 ) 

= 

OdO 

topz ( 1 ) 

= 

OdO 

topz (2 ) 

= 

-IdO 

topz (3) 

= 

-nh* ( (omz*ra/r0) **alfa) 

topz ( 4 ) 

= 

-IdO 

topz  ( 5 ) 

= 

(rhl/wj /rz/rz  -  ldO/rz  -  IdO /rhl ) *dexp ( - 

rhl / rz ) 

topz  ( 6 ) 

= 

(rh2/wj /rz/rz  -  ldO/rz  -  IdO / rh2 ) *dexp ( - 

rh2 / rz ) 

topz  ( 7 ) 

= 

IdO 

c 


Compute 


pressure  and 
PP  ( i )  = 

dpdr  = 
dpdz  = 
do  nn  = 
PP  (i) 
dpdr 
dpdz 
enddo 
PP  ( i )  = 

dpdr  = 
dpdz  = 
dpde  = 


its  derivatives 
OdO 
OdO 
OdO 
1,7 

=  pp(i)  +  c (nn) *top (nn) 

=  dpdr  +  c (nn) * (bot*topr (nn) 
=  dpdz  +  c (nn) * (bot*topz (nn) 

pp ( i ) /bot 
dpdr/bot2 
dpdz/bot2 
IdO /bot 


botr*top (nn) ) 
botz*top (nn) ) 


c  Gas  phase  limit 
else 


c  Compute  pressure  and  its  derivatives 
pp(i)  =  wj*ra*ea 

&  +  aj*(ldO  -  wj *ra/rhl ) *dexp ( -rhl /ra) 
&  +  bj*(ldO  -  wj *ra/rh2 ) *dexp ( -rh2 /ra) 
&  +  wj*ra*(qdetO  +  eO) 


dpdr  =  wj*ea 

&  +  aj*(rhl/ra2  -  w j /ra  -  wj /rhl ) *dexp ( -rhl /ra) 
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& 

& 


& 

& 


+  bj*(rh2/ra2  -  w j /ra  -  wj /rh2 ) *dexp ( -rh2/ra) 
+  w j  * (qdetO  +  eO ) 

dpdz  =  aj*(rhl/ra  -  wj  -  wj *ra/rhl ) *dexp ( -rhl/ra) 

+  bj*(rh2/ra  -  wj  -  wj *ra/rh2 ) *dexp ( -rh2 /ra) 

+  ra*wj* (qdetO  +  eO) 

dpde  =  wj*ra 
endif 


else  if  (ieos  .eq.  3)  then 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Hayes-II/JWL  EOS 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


ra 

=  rp ( i ) 

ra2 

=  ra*ra 

ea 

=  eip ( i 

za 

=  zp(i) 

rz 

=  ra*za 

omz 

=  ldO  - 

c  Solid  phase  limit 

if  (za  .le.  ztoll)  then 


c 


Compute  pressure  and  its  derivatives 

pp(i)  =  gh* (ea  -  beta*rO/ra  -  t4* ( (ra/rO) **alfa) 
&  +  t7 ) 

&  +  hi /nh* (( (ra/rO ) **nh)  -  ldO) 


& 

& 


dpdr  =  beta*r0*gh/ra2  -  alfa*gh*t4* (ra** (alfa-ldO) ) 
/ (rO**alfa) 

+  hi /rO* ( (ra/rO )* *nhml ) 


& 

& 


dpdz  =  gh*ea  -  beta*rO *gh/ra  +  alfa*gh*t4 
*  ( (ra/rO ) **alf a) 

+  hl/nh*(ldO  -  nhpl* ( (ra/rO ) **nh) ) 
dpde  =  gh 


c  Mixed  phases 

else  if  (ztoll  .It.  za  .and.  za  .It.  ztol2)  then 


c  Evaluate  denominator  functions 

bot  =  omz/gh  +  ldO/wj/ra 
if  (bot  .It.  Id-10)  then 
write  (  * , * )  '  ' 

write (*,*)  '  Zero  denonimator  term.' 

write  (  * , * )  '  ' 

stop 
endif 

bot2  =  bot*bot 
botr  =  -Id0/wj/ra2 
botz  =  -IdO/gh 

c  Evaluate  numerator  functions 
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top  (1) 
top  (2) 
top ( 3 ) 
top (4) 
top ( 5 ) 
top ( 6 ) 
top  (7) 
top  ( 8 ) 

c  Compute  derivatives 
topr  ( 1 ) 
topr (2 ) 
topr ( 3 ) 
topr ( 4 ) 
topr  ( 5 ) 
topr ( 6 ) 
topr  ( 7 ) 
topr ( 8 ) 


rhl / rz ) 
rh2/ rz) 


topz  ( 1 ) 
topz  (2 ) 
topz  ( 3 ) 
topz  ( 4 ) 
topz  (5) 

topz ( 6 ) 

topz  ( 7 ) 
topz  ( 8 ) 


c  Compute  pressure  and 
PP  ( i )  = 
dpdr  = 
dpdz  = 
do  nn  = 
PP  (i) 
dpdr 
dpdz 
enddo 
PP  ( i )  = 
dpdr  = 
dpdz  = 
dpde  = 


c  Gas  phase  limit 
else 


=  ea 

=  omz  -  rO/ra 

=  (omz**nh) * ( (ra/rO ) **alfa) 

=  omz 

=  (ldO/wj/ra  -  za/rhl ) *dexp ( -rhl /rz ) 

=  (ldO/wj/ra  -  za/rh2 ) *dexp ( -rh2 /rz ) 

=  za 

=  (omz**nhpl )*( (ra/rO )* *nh)  +  za  -  ldO 

for  numerator  functions 
=  OdO 
=  r0/ra2 

=  alfa/rO* (omz**nh) * ( (ra/rO) ** (alfa-ldO) ) 

=  OdO 

=  (rhl/wj/rz  -  ldO/wj  -  ldO ) *dexp ( -rhl /rz ) /ra2 
=  (rh2/wj/rz  -  ldO/wj  -  ldO ) *dexp ( -rh2/rz ) /ra2 
=  OdO 

=  nh/rO* (omz**nhpl) * ( (ra/rO) **nhml) 

=  OdO 
=  -ldO 

=  -nh* ( (omz*ra/rO ) **alf a) 

=  -ldO 

=  (rhl/wj /rz/rz  -  ldO/rz  -  ldO /rhl ) *dexp ( - 
=  (rh2/wj /rz/rz  -  ldO/rz  -  ldO / rh2 ) *dexp ( - 
=  ldO 

=  ldO  -  nhpl* ( (ra/rO*omz) **nh) 

its  derivatives 
OdO 
OdO 
OdO 
1,8 

=  pp(i)  +  c (nn) *top (nn) 

=  dpdr  +  c (nn) * (bot*topr (nn)  -  botr*top (nn) ) 

=  dpdz  +  c (nn) * (bot*topz (nn)  -  botz*top (nn) ) 

pp ( i ) /bot 
dpdr/bot2 
dpdz/bot2 
ldO/bot 


c  Compute  pressure  and  its  derivatives 
pp(i)  =  wj*ra*ea 

&  +  aj*(ldO  -  wj *ra/rhl ) *dexp ( -rhl/ra) 

&  +  bj*(ldO  -  wj *ra/rh2 ) *dexp ( -rh2/ra) 

&  +  wj*ra*(qdetO  +  eO) 

dpdr  =  wj*ea 

&  +  aj*(rhl/ra2  -  w j /ra  -  wj /rhl ) *dexp ( -rhl /ra) 

&  +  bj*(rh2/ra2  -  w j /ra  -  wj /rh2 ) *dexp ( -rh2 /ra) 

&  +  wj*(qdetO  +  eO) 
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& 

& 


dpdz  =  aj*(rhl/ra  -  wj  -  wj *ra/rhl ) *dexp ( -rhl /ra) 
+  bj*(rh2/ra  -  wj  -  wj *ra/rh2 ) *dexp ( -rh2/ra) 
+  ra*wj*(qdetO  +  eO) 

dpde  =  wj*ra 

endif 

else 

write  ( * , * )  '  ' 

write  (*,*)  '  Unknown  EOS' 

write  (  * , * )  '  ' 

stop 
endif 


enddo 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Estimate  mixture  temperature  Hayes-II/JWL  EOS  only 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
item  =  0 

if  (ieos  .eq.  3)  then 
item  =  1 
dtkmx  =  OdO 
denmx  =  OdO 


c  First 


c 


c 


temperature  estimate 
do  i  =  1 , imax-1 

if  (zp(i)  .gt.  ztol2)  then 
omz  =  OdO 
else 

omz  =  ldO  -  zp ( i ) 
endif 

denm  =  cvs*omz  +  cvg*zp(i) 
if  (denm  .gt.  denmx)  denmx  =  denm 

del  =  OdO 
de2  =  OdO 
de3  =  OdO 
de4  =  OdO 
de5  =  OdO 
de6  =  OdO 

if  (zp(i)  .It.  0.999d0)  then 
if  (zp(i)  .It.  ztol2)  then 
rs  =  omz*rp ( i ) 

del  =  t4* ( ( (rs/rO) **alfa)  -  ldO) 
de2  =  beta* (ldO  -  rO/rs) 
endif 

if  (zp(i)  .gt.  O.OOldO)  then 
if  (zp(i)  .gt.  ztoll)  then 
rg  =  zp ( i ) *rp ( i ) 
de3  =  aj /rhl*dexp (-rhl/rg) 
de4  =  bj /rh2*dexp (-rh2/rg) 
de5  =  a j /rhl *dexp ( -rhl /rO ) 
de6  =  b j /rh2 *dexp ( -rh2 /rO ) 
endif 


& 

& 


numr  =  eip(i)  -  omz* (del  -  de2) 

-  zp(i)*(de3  +  de4  -  de5  -  de6  -  qdetO 
+  eOcr) 


dtk(i)  = 

numr/ denm 

c 

write ( * , * 

)  '  zp  = 

' , zp  (i 

c 

write ( * , * 

)  '  del  = 

'  ,  del 

c 

write ( * , * 

)  '  de2  = 

'  ,  de2 

c 

write ( * , * 

)  '  de3  = 

'  ,de3 

c 

write ( * , * 

)  '  de4  = 

'  ,  de4 

c 

write ( * , * 

)  '  de5  = 

'  ,  de5 
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c 

c 


'  , dtk ( i ) 


write  (*,*)  '  de6  =  ' ,de6 

write  (*,*)  '  numr  =  ' ,numr,  1  dtk  = 

c  pause 

if  (dtk(i)  .It.  dtkmx)  dtkmx  =  dtk(i) 
enddo 

c  Check  the  temperature  difference  (is  T  <  TO?) 
if  (dtkmx  .It.  OdO)  then 
item  =  -1 

c  Calculate  the  internal  energy  correction  (fwded  to  next  time  level) 
eOcr  =  dtkmx*denmx/eta 

c  Apply  the  temperature  correction 
do  i  =  1 , imax-1 

omz  =  ldO  -  zp ( i ) 
denm  =  cvs*omz  +  cvg*zp(i) 
dtk(i)  =  dtk(i)  -  eOcr/denm 
enddo 
endif 

c  Calculate  the  corrected  temperature  field 
do  i  =  1 , imax-1 

tk(i)  =  tkO  +  dtk(i)  -  dtk(imax-l) 
c  tk(i)  =  dtk(i) 

enddo 

endif 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Update  particle  properties  and  positions 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
if  (ipar  .eq.  1)  then 

do  np  =  l,npar 

c  Compute  Reynolds  number 

ra  =  rp (peel (np) )* zp (peel (np) ) 
delu  =  up (peel  (np))  -  pu (np) 

adelu  =  dabs (delu) 

if  (adelu  .It.  Id-10)  adelu  =  ld-10 
rep  =  dip*ra*adelu/mu 


c 

write  (  * , * )  ' 

rep  =  ' , 

rep 

c 

if  (rep  .le. 

OdO)  then 

c 

write  (  * , * ) 

I  1 

c 

write ( * , * ) 

'  Rep  <= 

0  !  ' 

c 

write ( * , * ) 

'  cell  = 

' , peel (np) 

c 

write  (  * , * ) 

'  rp  =  ' 

, rp (peel (np) 

c 

write ( * , * ) 

'  zp  =  ' 

, zp (peel (np) 

c 

write ( * , * ) 

'  ra  =  ' 

,  ra 

c 

write ( * , * ) 

'  delu  = 

' , delu 

c 

write ( * , * ) 

'  adelu 

=  ' , adelu 

c 

write  (  * , * ) 

'  rep  = 

'  ,  rep 
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c 

c 


write  (  * , * )  '  ' 

stop 

c  endif 

c  Compute  particle  accelerations 
if  (idrg  .eq.  0)  then 

c  Spray  drag  law 

if  (rep  .It.  Id-10)  then 
cdp  =  OdO 

else  if  (rep  .le.  Id3)  then 

cdp  =  24d0/rep* (ldO  +  (rep**c23) /6d0) 
else 

cdp  =  0.44d0 
endif 

pa (np)  =  c316*mu*cdp*rep/rop/rdp/rdp*delu 

else  if  (idrg  .eq.  1)  then 

c  Rocket  drag  law 

if  (rep  .It.  Id-10)  then 
cdl  =  OdO 
cd2  =  OdO 
else 

cdl  =  24d0/rep  +  4 . 4d0/dsqrt (rep)  +  0.42d0 
cd2  =  c43* (1 . 75d0  +  150d0*alf21/rep) /alfl 
endif 

if  (alf2  .le.  0.08d0)  then 
cdO  =  cdl 

else  if  (0.08d0  .It.  alf2  .and.  alf2  .It.  0.45d0) 
cdO  =  (0 . 45d0-alf2) *cdl  +  (alf2-0 . 08d0) *cd2 
cdO  =  cd0/0.37d0 
else  if  (alf2  .gt.  0.45d0)  then 
cdO  =  cd2 
endif 

c  Mach  correction 

if  (imach  .eq.  1)  then 

mach  =  (adelu/ap (i) ) **4 . 63d0 
cdp  =  cd0*(ld0  +  dexp (-0 . 427d0/mach) ) 
else 

cdp  =  cdO 
endif 

pa (np)  =  cl8*pi*dip*dip*cdp*ra*adelu*delu/p0mas 
else 

write  ( * , * )  '  ' 

write  (*,*)  '  Unknown  drag  law.' 

write  (  * , * )  '  ' 

stop 


endif 

c 

write ( * , * ) 

'  rep  = 

'  ,  rep 

c 

write ( * , * ) 

'  cdp  = 

'  ,  cdp 

c 

write ( * , * ) 

'  pa  = 

' , pa (np) 
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then 


c 


write  (  * , * )  '  ' 


c  Compute  particle  velocity 

pup (np)  =  pu (np) 

c  Compute  particle  position 

pxp(np)  =  px (np) 


+  dt*pa (np) 

+  dt*pup(np) 


c  Set  default  particle  temperature 
ptkp(np)  =  tkO 


c  Update  particle  heat  transfer  and  temperature 
if  (ieos  .eq.  3)  then 


c  Compute  the  Nusselt  number  based  on  particle  Reynolds  number 
if  (rep  .le.  2d2)  then 

nup  =  2d0  +  0 . 106d0*rep*crppr 
else 

nup  =  2.274d0  +  0 . 6d0* (rep**0 . 76d0 ) *crppr 
endif 


c  Compute  the  heat  transfer  coefficient 
hp  =  tcon*nup/dip 

c  Compute  the  heat  transfer  coupling  term 

pq(np)  =  hp*pi*dip*dip* (tk (peel (np) )  -  ptk(np)) 

c  Compute  the  particle  temperature  change 
dtp  =  dt*pq(np) 

ptkp(np)  =  ptk(np)  +  dtp 


endif 


c  Check  particle  bounds 


c 

if  (pxp(np) 

.It.  xl)  then 

c 

pxp(np)  = 

xl 

c 

write  (  * , * ) 

i  '  Particle  ' , np. 

c 

stop 

c 

endif 

c 

if  (pxp(np) 

. gt.  x2)  then 

c 

pxp(np)  = 

x2 

c 

write  (  * , * ) 

i  '  Particle  ' , np. 

c 

stop 

c 

endif 

out  of  bounds . ' 


out  of  bounds . ' 


enddo 

endif 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Update  time  and  iteration  number 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
n  =  n  +  1 
time  =  time  +  dt 

write (*,*)  nstart+n, '  ' ,dt, '  ' ,time, '  ' , item 

write  (*,*)  ' pum  =  ' , pum 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  Solution  and  restart  file  output 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
if  (mod(n,ndmp)  .eq.  0)  then 
nfil  =  nfil  +  1 

c  Solution  file 

90  f ormat (' sol_' , i3 . 3, '. data ' ) 
write ( f ilex, 90 )  nfil 

open (22 , f ile=f ilex, f orm= ' formatted ' ) 
write (22,*)  '#  ' , time 

do  i  =  1 , imax-1 

xc  =  cl2*(x(i)  +  x(i+l)) 

write  (22,72)  xc,  rp  (i)  ,  up  (i)  ,  pp  (i)  ,  zp  (i)  ,  eip  (i)  ,  ap  (i)  , 
&  rxr ( i ) , tk ( i ) 

enddo 
close  (22 ) 

c  Particle  file 

91  f ormat (' par_ ', i3 . 3 data ' ) 
if  (ipar  .eq.  1)  then 

write (parex, 91 )  nfil 

open (22 , f ile=parex, form= ' formatted ' ) 
do  np  =  l,npar 

write (22,*)  pxp(np),'  ' ,pup(np),'  ' ,ptkp(np) 
enddo 
close  (22 ) 
endif 

c  Derivatives  file 

open (22 , f ile= ' deriv . data ' , f orm= ' formatted ' ) 
do  i  =  1 , imax-1 

write(22,*)  i,'  ' , derv ( i , 1 ) , '  ',derv(i,2) 

enddo 
close  (22 ) 

c  L/R  Z  files 

c  open (22 , f ile= ' zlzr . data ' , form= ' formatted ' ) 

c  do  i  =  1 , imax 

c  write(22,*)  i,'  ',zzl(i),'  ',zzr(i) 

c  enddo 

c  close  (22 ) 


c  Restart  file 

open ( 40 , f ile= ' restart . data ' , f orm= ' unformatted ' ) 

write (40)  nstart+n 

write (40)  nfil 

write (40)  time 

do  i  =  1 , imax-1 

write  (40)  rp ( i ) , pp ( i ) , up ( i ) , zp ( i ) 
enddo 
close (40) 

endif 

c  Reset  arrays 
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i  i  = 

1 , imax-1 

r  (i) 

=  rp ( i ) 

u(i) 

=  up ( i ) 

z  (i) 

=  zp(i) 

ei  (i) 

=  eip ( i 

P  ( i ) 

=  PP ( i ) 

a  (i) 

=  ap(i) 

enddo 

92  format (2x,dl5.9,2x,dl5.9,2x,dl5.9,2x,dl5.9,2x,i5) 

if  (ipar  .eq.  1)  then 
do  np  =  l,npar 
px(np)  =  pxp(np) 
pu (np)  =  pup (np) 
ptk(np)  =  ptkp(np) 

if  (idbgp  .eq.  1)  write (110+np, 92)  time, pxp (np) , 

&  pup (np) , pa (np) , peel (np) 

enddo 

endif 

c  pause 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c  End  of  solver  loop 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

enddo 

c  Termination  codes 

if  (  time  .gt.  tend)  then 
write ( * , * )  '  ' 

write (*,*)  '  TIME  >  TEND.' 

else  if  (n  .ge.  nstp)  then 
write ( * , * )  '  ' 

write  (*, *)  '  N  >  NSTP.  ' 

else 

write (*,*)  '  UNKNOWN  TERMINATION  CRITERIA.' 

endif 

c  End  of  main  program 
stop 
end 
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